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Abstract 

A  new  theory  and  algorithm  for  scatterer  classification  in  synthetic  aperture  radar 
(SAR)  imagery  is  presented.  The  automated  classification  process  is  operationally 
efficient  compared  to  existing  image  segmentation  methods  requiring  human  super¬ 
vision. 

The  algorithm  reconstructs  coarse  resolution  subimages  from  subdomains  of  the 
SAR  phase  history.  It  analyzes  local  peaks  in  the  subimages  to  determine  locations 
and  geometric  shapes  of  scatterers  in  the  scene.  Scatterer  locations  are  indicated 
by  the  presence  of  a  stable  peak  in  all  subimages  for  a  given  subaperture,  while 
scatterer  shapes  are  indicated  by  changes  in  pixel  intensity.  A  new  multi-peak  model 
is  developed  from  physical  models  of  electromagnetic  scattering  to  predict  how  pixel 
intensities  behave  for  different  scatterer  shapes.  The  algorithm  uses  a  least  squares 
classifier  to  match  observed  pixel  behavior  to  the  model.  Classification  accuracy 
improves  with  increasing  fractional  bandwidth  and  is  subject  to  the  high-frequency 
and  wide-aperture  approximations  of  the  multi-peak  model. 

For  superior  computational  efficiency,  an  integrated  fast  SAR  imaging  technique  is 
developed  to  combine  the  coarse  resolution  subimages  into  a  final  SAR  image  having 
fine  resolution.  Finally,  classification  results  are  overlaid  on  the  SAR  image  so  that 
analysts  can  deduce  the  significance  of  the  scatterer  shape  information  within  the 
image  context. 
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PHASE  HISTORY  DECOMPOSITION  FOR  EFFICIENT  SCATTERER 


CLASSIFICATION  IN  SAR  IMAGERY 

I.  Introduction 

This  chapter  provides  an  introduction  and  overview  of  the  dissertation  document. 
Section  1.1  is  an  executive  summary  of  the  research  motivation,  hypotheses,  and 
findings.  Section  1.2  highlights  the  operational  needs  which  are  addressed  by  this 
dissertation.  Section  1.3  lists  the  conference  papers,  journal  articles,  and  other  deliv¬ 
erables  which  have  been  produced  in  connection  with  this  dissertation.  Section  1.4 
wraps  up  with  an  outline  of  the  document  and  the  organization  of  its  chapters. 

1.1  Executive  Summary 

Resource  management  is  an  ongoing  need  in  defense  operations.  As  a  result, 
synthetic  aperture  radar  (SAR)  imaging  and  classification  algorithms  are  needed  to 
rapidly  queue  human  operators  and  precision  algorithms  to  regions  of  high  inter¬ 
est.  This  dissertation  describes  a  new  SAR  imaging  and  classification  theory  as  a 
foundation  from  which  to  build  rapid  queuing  solutions  for  improved  operational  ef¬ 
ficiency.  The  theory  is  demonstrated  in  a  new  algorithm  based  on  efficient  imaging 
and  classification  techniques. 

Phase  history  decomposition  is  a  highly  efficient  technique  for  SAR  image  re¬ 
construction,  where  subimages  are  produced  as  an  intermediate  step  of  the  imaging 
process.  The  subimages  have  coarser  resolution  than  the  final  image,  but  have  been 
shown  to  provide  insight  into  the  anisotropic  and  dispersive  nature  of  objects  in  the 
image  scene.  While  research  on  the  anisotropic  nature  of  scatterers  in  subimages 
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has  been  on-going  for  some  time,  research  on  the  dispersive  nature  of  scatterers  in 
subimages  is  quite  nascent  and  worthy  of  additional  investigation.  In  response,  this 
dissertation  investigates  the  concept  of  scatterer  classification  by  subimage  analysis 
and  develops  a  fundamental  and  holistic  theory  for  this  emerging  research  area.  This 
research  hypothesizes  that 

•  it  is  possible  to  locate  and  classify  canonical  scatterers  by  observing  the  inten¬ 
sities  of  subimage  pixels,  and 

•  phase  history  decomposition  makes  this  approach  to  classification  highly  effi¬ 
cient. 

The  key  findings  of  this  research  effort  include 

1.  image  peaks  due  to  a  distributed  canonical  scatterer  can  be  modeled  with  a 
simple  equivalent  canonical  point  scatterer  [58], 

2.  the  intensities  of  subimage  peaks  reveal  the  locations  and  types  of  canonical 
point  scatterers  in  a  SAR  scene  [57], 

3.  the  classification  and  imaging  errors  associated  with  phase  history  decomposi¬ 
tion  are  controllable  [58,  59],  and 

4.  the  proposed  approach  is  novel,  efficient,  and  foundational  [58,  59]. 

The  first  key  finding  results  from  development  and  study  of  a  new  scattering  model 
called  the  multi-peak  model.  The  second  key  finding  results  from  a  new  scatterer 
classification  algorithm  called  the  spectrum  parted  linked  image  test  (SPLIT).  Imag¬ 
ing  accuracy  in  the  third  key  finding  and  computational  efficiency  in  the  fourth  result 
from  a  new  integrated  algorithm  that  combines  fast  SAR  imaging  techniques  with 
scatterer  classification.  The  greatest  benefit  of  the  new  theory  is  the  operational  ef¬ 
ficiency  derived  by  automatically  displaying  scatterer  classification  results  within  the 
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context  of  the  SAR  image.  This  is  a  significant  improvement  over  existing  scatterer 
classification  methods  [3,  30,  46,  47,  48,  62,  73]  that  require  human  supervision  to 
ensure  accurate  classification. 

1.2  Operational  Needs 

One  of  the  ongoing  needs  in  defense  operations  is  resource  management  due  to  the 
fact  that  there  exists  more  imagery  collection  capability  than  processing  capability 
[13,  19,  32,  102,  115].  Therefore,  SAR  imaging  algorithms  are  needed  to  provide  rapid 
queuing  of  operators  and  precision  algorithms  to  regions  of  high  interest  [63,  142],  In 
general,  precision  algorithms  are  not  adequate  in  meeting  this  need  because  they  are 
computationally  inefficient,  operationally  inflexible,  or  both  [54,  80,  119,  127]. 

For  instance,  precision  SAR  imaging  algorithms  improve  SAR  image  quality  by 
using  better  geometric  approximations  commonly  made  in  the  imaging  algorithm 
[24,  36,  76].  However,  increasing  the  order  of  these  geometric  approximations  comes 
at  an  increased  computational  cost  [24,  36].  Therefore,  computational  resource  man¬ 
agement  is  optimized  when  use  of  precision  SAR  imaging  algorithms  is  limited  to 
regions  where  increased  precision  is  warranted. 

Precision  target  recognition  algorithms  are  notoriously  sensitive  to  operational 
conditions,  which  cause  them  to  be  inflexible  outside  of  a  specific  operational  scenario 
[54,  80,  119,  127].  In  this  case,  the  precision  algorithms  must  be  used  selectively  and 
queued  by  experienced  analysts  based  on  operational  parameters  and  image  context. 
Unfortunately,  this  creates  an  operational  bottleneck  by  demanding  human  resources 
be  used  to  manage  precision  target  recognition  algorithms.  This  is,  in  effect,  the  exact 
opposite  of  what  is  needed  to  improve  resource  management  in  defense  operations. 

In  order  to  provide  rapid  queuing  of  operational  resources,  it  is  acceptable  for  a 
SAR  imaging  and  classification  algorithm  to  sacrifice  some  precision  in  order  to  obtain 


3 


*  ✓ 
\ 


\  ' 

. . 

Integrated  Scatterer 
Classification  and  Imaging  by 
Phase  History  Decomposition 


Figure  1.  The  dissertation  research  improves  theoretical  knowledge  in  two  new  areas. 


efficiency  and  flexibility.  In  keeping  with  this  principle,  this  dissertation  describes  a 
new  SAR  imaging  and  classification  theory  as  a  foundation  from  which  to  build  rapid 
queuing  solutions  for  improved  operational  efficiency. 

1.3  Contributions 

This  dissertation  research  improves  theoretical  knowledge  in  two  new  areas,  as 
shown  in  Figure  1.  It  develops  new  theory  for  Scatterer  Classification  by  Phase 
History  Decomposition  and  combines  this  with  existing  theories  in  Scattering  Matrix 
Decomposition  and  Domain  Decomposition  Imaging.  It  also  develops  unique  design 
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principles  to  address  the  particular  difficulties  associated  with  combining  these  into 
an  integrated  theory  and  algorithm. 

The  research  presented  in  this  dissertation  resulted  in  three  published  conference 
papers  and  two  journal  article  submissions  [55,  56,  57,  58,  59].  In  addition,  a  fast  con¬ 
volution  backprojection  and  scatterer  classification  code  was  fully  developed  and  will 
be  provided  to  offices  who  sponsored  elements  of  this  work.  It  will  also  be  considered 
for  integration  into  the  AFIT  LORE  Processing  INtegrated  Environment  (ALPINE). 
From  this  code,  alternate  versions  of  the  polar  format  algorithm,  range  Doppler  al¬ 
gorithm,  and  omega-k  algorithm  were  also  developed  for  experimental  purposes. 

Last,  this  dissertation  produces  a  theoretical  foundation  for  follow-on  research  in 
the  following  areas: 

•  discontinuous  phase  histories, 

•  improved  parameter  estimation  using  advanced  detection  and  estimation  theory, 

•  extension  to  bi-static  and  3D  SAR, 

•  blended  domain  decomposition  and  decimation  techniques,  and 

•  additional  uses  for  the  multi-peak  model  and  SPLIT  algorithm. 

These  are  described  in  more  detail  in  Subsection  7.2  of  the  Conclusion. 

1.4  Organization 

The  dissertation  is  organized  as  follows.  Chapter  II  presents  a  survey  of  the  ma¬ 
ture  research  areas  in  Figure  1.  Using  a  combination  of  tutorial  and  literature  review, 
it  presents  theory  and  trends  in  the  topics  of  SAR  imaging,  domain  decomposition, 
canonical  scattering  models,  and  scatterer  classification.  Chapter  III  presents  the 
research  objective  in  general  terms.  It  serves  as  a  transition  between  the  background 
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section  and  the  more  detailed  theoretical  developments  and  research  findings  con¬ 
tained  in  the  remaining  sections.  Chapter  IV  develops  the  new  multi-peak  model 
for  canonical  scatterers.  It  describes  how,  under  a  wide-angle  condition,  the  imaging 
process  integrates  out  the  azimuth  dependency  of  distributed  canonical  scatterers.  In 
this  way,  image  peaks  due  to  a  distributed  canonical  scatterer  can  be  modeled  as  due 
to  an  equivalent  canonical  point  scatterer.  Chapter  V  develops  the  SPLIT  algorithm, 
which  uses  subimage  pixel  intensities  to  estimate  the  locations  of  canonical  scatter¬ 
ers  as  well  as  their  frequency  dependencies.  SPLIT  classifies  the  observed  canonical 
scatterers  using  the  multi-point  model  to  deduce  the  likelihood  that  a  certain  type  of 
scatterer  is  present.  Chapter  VI  develops  the  integrated  algorithm  which  combines 
SPLIT-based  classification  with  domain  decomposition  imaging.  The  integrated  al¬ 
gorithm  is  shown  to  be  efficient  in  that  it  provides  scatterer  classification  information 
without  increasing  the  computational  complexity  of  SAR  imaging  algorithms.  The 
combined  results  provide  more  information  about  the  scene  than  a  SAR  image  can 
provide  alone.  Chapter  VII  concludes  with  an  overview  of  the  key  hirelings  and  con¬ 
tributions  in  this  dissertation,  a  summary  of  the  advantages  and  limitations  of  the 
algorithm,  and  suggestions  for  follow-on  research. 
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II.  Background 


This  chapter  provides  the  background  needed  to  understand  the  theory  and  trends 
in  areas  related  to  this  research.  It  serves  as  a  tutorial  and  includes  reviews  of  the 
seminal  and  current  literature,  where  appropriate,  along  with  observations  of  trends. 
Section  2.1  begins  with  an  overview  of  basic  SAR  imaging  concepts.  It  builds  on  these 
basic  concepts  to  present  the  more  advanced  concepts  of  phase  history  decomposition 
imaging  and  the  sum  of  scattering  centers  model.  Section  2.2  presents  the  set  of 
canonical  scatterers  used  to  obtain  a  parsimonious  sum  of  scattering  centers.  It 
explains  how  the  amplitude  and  polarization  responses  of  canonical  scatterers  are 
derived  from  physical  models  of  electromagnetic  scattering.  Section  2.3  gives  an 
overview  of  some  basic  principles  of  feature  extraction  and  classification,  to  include 
the  least  squares  classifier  featured  in  this  dissertation.  It  concludes  with  an  example 
of  least  squares  classification  using  the  polarimetric  model  parameters  for  canonical 
scatterers. 

2.1  SAR  Phase  History 

SAR  images  are  typically  reconstructed  from  a  SAR  phase  history,  which  repre¬ 
sents  the  SAR  signal  in  the  spectral  domain.  This  section  presents  how  the  SAR 
phase  history  is  originated.  Then  it  presents  principles  of  SAR  imaging,  including 
domain  decomposition  imaging.  Finally,  it  presents  how  the  phase  history  can  be 
modeled  as  a  sum  of  phase  histories  clue  to  multiple  scattering  centers. 

2.1.1  Origination. 

Monostatic  SAR  systems  measure  the  electromagnetic  reflectivity  of  objects  in 
the  radar  held  of  view  and  rely  on  subsequent  signal  processing  to  reconstruct  an 
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Surveillance  Area 


Figure  2.  Notional  airborne  SAR  system. 
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Figure  3.  Notional  phase  history  domain. 


estimate  of  the  reflectivity.  The  estimate  is  presented  as  a  SAR  image  recovered 
from  quadrature-demodulated  samples  of  the  backscattered  electric  field  received  at 
discrete  and  different  frequencies  and  aspect  angles  along  the  synthetic  aperture  [24, 
36,  37,  76,  133].  These  discrete  samples  are  collectively  referred  to  as  the  SAR  phase 
history.  The  term  ‘phase  history’  refers  to  the  phase  differences  corresponding  to  the 
relative  locations  of  each  scatterer  in  the  scene. 

A  notional  airborne  SAR  system  is  depicted  in  Figure  2  where  discrete  samples 
of  the  scattered  electric  field  are  collected  over  the  flight  path.  The  phase  history 
is  typically  displayed  as  a  manifold  at  sample  coordinates  in  the  spectral  domain 
[24,  76].  For  example,  a  notional  SAR  phase  history  domain,  which  states  G  is 
a  function  of  (/,  9),  is  depicted  in  Figure  3,  where  the  ~  symbol  denotes  that  the 
phase  history  is  a  complex-valued  function.  The  discrete  samples  are  contained  in 
the  domain  f2  =  [f,  0],  where  f  is  a  vector  of  sample  points  in  frequency  and  6  is 
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a  vector  of  sample  points  in  azimuth.  For  convenience,  the  individual  phase  history 
samples  are  not  shown.  This  allows  drawings  of  the  SAR  phase  history  to  be  readily 
differentiated  from  drawings  of  SAR  imagery  which  explicitly  show  individual  pixels. 
This  convention  is  used  throughout  the  dissertation. 

The  sample  frequencies  are  determined  by  the  transmitted  electromagnetic  field, 
usually  produced  with  a  linear  frequency  modulated  (LFM)  signal  [24,  76].  For  mul¬ 
tiple  polarization  channels,  a  separate  phase  history  is  produced  for  each  channel.  A 
particularly  important  aspect  of  the  SAR  phase  history  is  that  it  has  frequency  and 
azimuth  diversity,  and  in  the  case  of  polarimetric  SAR,  it  has  polarimetric  diversity  as 
well.  The  remaining  discussion  assumes  that  the  sampling  rate  is  sufficient  to  prevent 
aliasing  in  the  image  and  that  all  amplitude  variations  due  to  antenna  gain  pattern 
and  spherical  wave  propagation  are  normalized  between  samples. 

2.1.2  SAR  Imaging. 

The  image,  g,  is  reconstructed  using  an  appropriate  transformation  from  the  spec¬ 
tral  domain  to  the  spatial  domain.  Recall  that  the  phase  history,  G,  is  defined  over 
finite  regions  of  support,  where  fc  is  the  center  frequency  of  the  phase  history  with 
bandwidth  R,  and  9C  is  the  center  angle  of  the  phase  history  with  aperture  width 
0.  In  this  case,  the  finite  regions  of  support  can  be  represented  by  a  band-limited 
filter  or  window  in  frequency,  HP(f  —  /c),  and  an  aperture-limited  filter  or  window  in 
azimuth,  Hq(9  —  9c),  where  the  windows  have  region  of  support  Hs{f)  G  [— B/2,  B / 2] 
and  Hq(9)  G  [—0/2,  0/2],  respectively.  Hence,  the  2D  image  is  reconstructed  from 
the  windowed  phase  history  as 

g(x,y)  =  B  -  Jc)He($  -  0c)G(j,  0)}  ,  (1) 
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where  (x,y)  are  image  coordinates  and  £>{•}  is  a  Fourier-based  imaging  operator 
that  maps  from  (/,  6)  to  (x,y).  The  choice  of  Hb  and  Hq  depends  on  the  need  for 
image  resolution  versus  image  contrast.  For  convenience,  the  g  dependence  on  Hb 
and  Hq  is  suppressed  in  the  notation.  Although  the  phase  history  and  image  are 
digitally  sampled  and  processed  in  practice,  the  variables  /,  6 ,  x,  and  y  are  expressed 
as  continuous  for  ease  of  notation.  Note  that  if  multiple  polarization  channels  are 
available,  the  transformation  is  performed  for  each  channel’s  phase  history,  resulting 
in  a  set  of  usually  two  to  four  polarization  diverse  images. 

2. 1.2.1  Subimages. 

The  relationship  between  spectral  bandwidth  in  the  phase  history  and  spatial 
resolution  in  the  SAR  image  is  a  manifestation  of  the  Gabor  limit  [60].  In  short, 
the  spatial  resolution  is  inversely  proportional  to  the  spectral  bandwidth.  Therefore, 
it  is  possible  to  produce  coarse  resolution  subimages  from  subdomains  of  the  phase 
history.  An  example  of  this  is  illustrated  in  Figure  4  where  Figure  4(b)  shows  a  coarse 
resolution  subimage.  In  this  case,  the  subwindows,  Hb'  and  Hq>,  decompose  the 
spectral  domain,  where  the  regions  of  support  for  these  subwindows,  B'  <  B  and 
0'  <  0,  are  called  subbands  and  subapertures,  respectively.  Multiple  subdomains 
may  be  created  by  simply  shifting  the  subwindows  to  a  discrete  number  of  center 
frequencies.  In  this  case,  center  frequencies  are  annotated  by  subscripts  i  and  j,  and 
the  short  hand  notations  Hbh  =  Hs'(f  —  fCi)  and  Hq>3  =  Hq^IO  —  9Cj)  are  used 
throughout  the  dissertation,  where  convenient.  Thus,  the  reconstructed  subimages 
are  annotated  accordingly  as 

gij(x,y)  =  B{HB,iH&jG{f,e)}.  (2) 
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(a)  A  fine  resolution  SAR  image  is  reconstructed  from  the  full  domain  of  the 
phase  history. 


A  G(f  0)  Hb  -iHe  jGif,  0)  g(x,  y ) 


(b)  A  coarse  resolution  SAR  subimage  is  reconstructed  from  a  subdomain  of  the 
phase  history. 


Figure  4.  Phase  history  decomposition  produces  multiple  coarse  resolution  subimages. 


When  appropriate,  the  subscript,  p,  can  be  added  to  gl3  and  G  to  denote  subimages 
produced  from  phase  histories  obtained  from  different  polarization  channels. 


2. 1.2. 2  Phase  History  Decomposition. 

The  phase  history  can  be  replicated  by  using  a  series  of  subwindows  designed 
and  weighted  so  that  their  summation  approximates  the  desired  fullband  and  full 
aperture  windows.  The  summations  are  expressed  as  HB  ~  Y^iciHB'(f  —  fa)  and 
Hq  ~  cjHefd  —  6Cj),  where  ct  and  c3  are  the  weights.  The  shorthand  notations 
H&  =  Hq(9  —  6C )  and  HB  =  HB(f  —  fc )  are  used  here  and  throughout  the  disserta¬ 
tion,  as  appropriate.  In  this  case,  the  resulting  coarse  resolution  subimages  can  be 
interpolated  to  a  finer  resolution  and  summed,  where  the  result  approximates  the  fine 
resolution  image  conventionally  reconstructed  from  the  full  domain  of  the  phase  his- 
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Figure  5.  Examples  of  domain  decomposition  imagery  of  four  automobiles  in  a  parking 
lot  taken  from  the  Gotcha  public  release  data  set. 


tory  [12].  This  assumes  a  linear  imaging  operator,  as  is  common  [76],  and  is  expressed 
as 

9  =  ^c3^CiI{9ij}  ■  (3) 

3  * 

where  Z{-}  is  the  interpolation  operator  and  the  ^  symbol  indicates  that  the  result 
is  an  approximation  to  the  full  domain  image  in  Equation  (1).  This  process  is  called 
domain  decomposition  imaging,  where  typical  implementations  have  a  controllable 
error  in  the  approximation  [12].  Examples  of  conventional  imaging  and  domain  de¬ 
composition  imaging  using  cubic  interpolation  are  shown  in  Figures  5(a)  and  5(b), 
respectively.  In  this  case,  error  in  the  approximation  is  sufficiently  controlled  so  that 
the  images  are  visually  indistinguishable.  The  scene  consists  of  four  automobiles  in  a 
parking  lot  taken  from  the  Gotcha  public  release  data  set  [25]. 

Domain  decomposition  imaging  is  traditionally  motivated  by  the  desire  to  reduce 
the  computational  complexity  of  certain  imaging  algorithms  by  accepting  a  control¬ 
lable  error  in  imaging  accuracy  [12].  Efficiency  is  attained  by  an  iterative,  multi-level 
decomposition  and  aggregation  of  subimages.  The  overall  computational  complexity 
of  multi-level  domain  decomposition  algorithms  is  0(N2  log  N),  for  an  N  x  N  image 
[12].  Domain  decomposition  can  be  implemented  with  any  of  the  conventional  SAR 
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Table  1.  Computational  Complexity  of  Conventional  SAR  Imaging  Algorithms  [24,  36, 
76]. 


Mode 

Algorithm 

Complexity 

spotlight 

Matched  Filter 

0(A4) 

Convolution  Backprojection 

0(N 3) 

Polar  Reformating 

0(A2  log  A) 

stripmap 

Range  Doppler 

0(A2  log  A) 

Chirp  Scaling 

0(N2  log  A) 

Omega-K 

0(A2  log  A) 

imaging  algorithms  listed  in  Table  1.  However,  the  benefit  of  a  reduced  order  of 
computational  complexity  will  only  be  realized  for  the  matched  filter  and  convolution 
backprojection  algorithms. 

2. 1.2. 3  Observations. 

SAR  imaging  is  a  mature  area  of  research  with  many  textbooks  dedicated  to  the 
various  methods  and  their  applications.  While  the  matched  filter  (MF)  algorithm 
provides  the  most  flexibility  and  best  image  quality  of  the  algorithms,  its  computa¬ 
tional  complexity  of  0(N 4)  is  exceptionally  high  compared  to  other  methods.  The 
other  algorithms  obtain  computational  efficiency  by  use  of  batch  processes  with  geo¬ 
metric  approximations  which  limit  imaging  accuracy  and  flexibility.  However,  these 
limitations  can  be  managed  so  that  they  are  insignificant  for  most  SAR  applications. 
The  algorithms  with  computational  complexity  of  0(A2  log  A)  obtain  computational 
efficiency  through  2D  Fourier  transforms  of  rectangular  formatted  phase  history  data. 
In  contrast,  the  convolution  backprojection  (CBP)  algorithm  uses  the  projection  slice 
theorem  with  polar  formatted  phase  histories,  resulting  in  a  computational  complex¬ 
ity  of  0(A3)  with  higher  order  [40].  Because  of  this,  CBP  is  usually  employed  only 
when  its  superior  flexibility  in  choosing  the  locations  of  image  pixels  is  needed.  An 
example  situation  where  such  flexibility  is  desired  is  the  case  of  2D  imaging  of  the 
surface  of  the  earth  over  very  wide-angle  or  full  360°  apertures.  In  this  case,  the  CBP 
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algorithm  allows  for  image  pixel  locations  to  be  selected  to  match  the  digital  terrain 
elevation  data  (DTED)  of  the  scene  for  superior  image  quality. 

The  fast  SAR  imaging  techniques  ensure  a  computational  complexity  at  or  near 
0(N  log  N)  with  a  controllable  amount  of  image  artifacts  introduced  into  the  final 
image.  In  general,  the  fast  SAR  imaging  techniques  only  provide  a  computational 
advantage  for  the  CBP  and  MF  imaging  algorithms,  but  they  are  not  limited  to 
these.  There  are  two  primary  techniques  for  fast  SAR  imaging:  domain  decom¬ 
position  and  domain  decimation.  Domain  decomposition  imaging  produces  coarse 
resolution  subimages  with  diversity  in  frequency,  azimuth,  and  polarization,  which 
can  be  exploited  for  scatterer  classification.  As  such,  a  single  level  version  of  the 
multilevel  domain  decomposition  technique  in  Reference  [12]  is  used  throughout  this 
dissertation.  Alternately,  domain  decimation  produces  full  resolution  subimages  of 
limited  extent  that  are  diverse  in  location.  The  final  image  is  reassembled  from  these 
subimages  by  a  process  resembling  a  mosaic.  Because  the  pixel  locations  of  the  full 
resolution  subimages  can  be  adjusted  with  precision,  domain  decimation  is  the  pre¬ 
ferred  method  for  applications  where  the  image  pixels  are  matched  to  DTED,  such 
as  with  the  Gotcha  radar. 

It  is  conceivable  to  combine  the  decimation  and  decomposition  techniques,  al¬ 
though  this  has  not  been  reported  in  the  literature.  In  this  way,  multilevel  domain 
decimation  could  be  used  to  form  the  coarse  resolution  subimages  matching  DTED. 
Then,  the  coarse  resolution  subimages  can  be  used  for  scatterer  classification  and 
subsequent  single-level  domain  decomposition  can  be  used  to  form  the  final  image. 

2. 1.2. 4  Subaperture  Imaging. 

Subaperture  imaging  is  a  decomposition  of  the  phase  history  in  azimuth  only 
and  can  be  used  with  any  of  the  conventional  imaging  algorithms.  Its  use  is  typi- 
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cally  motivated  by  the  need  to  limit  azimuth  resolution,  reduce  speckle,  or  both.  In 
practice,  SAR  systems  are  often  less  restricted  in  angular  bandwidth  than  frequency 
bandwidth.  This  relationship  is  particularly  true  for  airborne  circular  SAR  systems, 
like  the  one  in  Figure  2  of  Section  2.1.1  which  can  orbit  a  given  scene  and  produce 
a  phase  history  sampled  over  a  full  360°  in  azimuth.  In  this  case,  subapertures  of 
the  data  are  processed  in  turn.  For  surveillance  applications,  a  series  of  subimages 
is  often  viewed  in  sequence  to  simulate  frames  of  streaming  video,  called  video  SAR 
[18,  109].  In  addition,  when  the  subimages  are  registered  to  a  common  grid,  they  can 
be  coherently  summed  to  produce  a  full  aperture  image  expressed  as 


g(x,y)  =  ^CjB  {HBH&jG(f,6)}}  =  cjgj{x}  y).  (4) 

j  j 

Figure  5(a)  is  an  example  of  subaperture  imaging.  Here,  multiple  subimages  are 
reconstructed  from  successive  2°  subapertures  and  are  summed  to  approximate  a  full 
360°  aperture  SAR  image. 

The  grainy  look  in  Figures  5(a)  and  5(b)  is  attributed  to  a  common  SAR  imaging 
effect  called  speckle.  Speckle  in  SAR  imagery  results  from  a  combination  of  having 
multiple  scatterers  and  only  finite  processing  resolution.  The  phase  histories  of  the 
unresolved  scatterers  produce  a  random  sum  with  characteristic  appearance  in  the 
SAR  image,  although  not  appearing  in  photographs  of  the  same  scene.  In  fact,  speckle 
has  been  shown  to  be  well-modeled  by  a  random  phase  process  [113].  By  taking  the 
root  mean  squared  (RMS)  of  pixel  values  over  multiple  subaperture  images,  areas 
with  highly  random  pixel  intensities  become  smoother  and  the  effects  of  speckle  are 
reduced.  This  is  also  called  multilook  imaging  and  is  expressed  as  [36] 


9{x,y) 


1 9j(x,y)f 


(5) 
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An  example  of  multilook  imaging  is  shown  in  Figure  5(c)  using  the  same  data  and 
subapertures  as  in  Figure  5(a).  Speckle-reduced  SAR  imagery  is  often  preferred  be¬ 
cause  it  can  be  easier  to  interpret  and  exploit,  even  for  the  case  when  the  subaperture 
images  have  coarser  resolution  than  the  full  aperture  image  [116].  Multilook  can  be 
combined  with  domain  decomposition  imaging  as 


9(x,y)  = 


\ 


£■ 


Ci^{dij} 


(6) 


where  the  ^  symbol  indicates  that  the  result  is  an  approximation  to  the  multilook 
image  in  Equation  (5).  An  example  of  multilook  domain  decomposition  imaging  using 
cubic  interpolation  is  given  in  Figure  5(d).  In  this  case,  error  in  the  approximation 
is  sufficiently  controlled  so  that  the  multilook  images  in  Figures  5(c)  and  5(d)  are 
visually  similar.  Note  that  the  imaging  process  is  repeated  for  each  polarization 
channel  separately.  The  resulting  set  of  polarization  diverse  images  can  be  non- 
coherently  summed,  if  desired. 


2. 1.2. 5  Observations. 

With  the  recent  availability  of  very-wide  angle  SAR  systems,  interest  has  increased 
in  exploiting  the  benefits  of  these  systems.  As  a  result,  many  different  methods  and 
techniques  for  combining  subimages  have  been  proposed.  However,  even  though  other 
non-coherent  summations  of  subaperture  images  have  been  discussed  in  the  literature 
[110],  the  domain  decomposition  Equations  (3)  and  (6)  are  commonly  used  in  practice, 
and  are  exclusively  used  throughout  this  dissertation. 
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Figure  6.  Vector  from  the  antenna  phase  center  to  the  qth.  scattering  center. 


2.1.3  Sum  of  Scattering  Centers  Model. 

By  linear  superposition,  a  SAR  phase  history  can  be  modeled  as  a  sum  of  phase 
histories  due  to  Q  scattering  centers  expressed  as  [138] 

Q- 1 

G(f,e)  =  Y,s,(f,e)e-W’K  (7) 

<7=0 


where  Sq(f,  9)  is  the  amplitude  function,  k  =  2nf/c  is  the  wavenumber  having  speed 
of  light,  c,  and  rq  is  a  vector  from  the  antenna  phase  center  to  the  qth  scattering 
center,  as  shown  in  Figure  6. 

A  linear  imaging  operator  is  commonly  used  in  practice,  and  in  this  case,  the  SAR 
image  can  be  modeled  as  a  linear  superposition  of  images  due  to  Q  scattering  centers 
expressed  as 

Q—i 

g(x,y)  =  ^2$q(x,y),  (8) 

g=0 


where 


sq(x,y)  =  B{HBHeSq(f,6)e-j2k^ '} 


(9) 


is  the  resultant  of  the  imaging  operator  acting  on  the  qth  windowed  phase  history  in 
Equation  (7). 
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Consider  a  scene  consisting  of  a  single  scattering  center.  In  this  case,  the  phase 
history  samples  contain  amplitude  and  phase  differences  corresponding  to  the  am¬ 
plitude  function  and  location  of  the  scatterer.  If  for  every  sample,  one  knows  the 
distance  from  the  SAR  system  to  a  reference  point  in  the  scene  (usually  the  scene 
center),  then  one  can  infer  the  location  of  the  scatterer  with  respect  to  the  refer¬ 
ence.  Furthermore,  by  applying  an  appropriate  phase  shift,  the  samples  will  integrate 
coherently. 

The  exact  locations  of  point  scatterers  are  usually  not  known  a  priori.  Therefore, 
the  imaging  operator  integrates  over  a  grid  of  locations,  each  corresponding  to  a 
unique  phase  shift.  This  grid  determines  the  pixel  locations  of  the  SAR  image,  and 
in  the  2D  case,  is  called  the  imaging  plane.  The  accuracy  of  the  coherent  integration, 
and  thus  the  accuracy  of  the  image,  is  limited  to  time-invariance  of  the  scatterers, 
field  geometry,  and  SAR  system  [76]. 

Typically  the  imaging  operator  assumes  a  scene  comprised  of  ideal  point  scatterers. 
In  this  case,  Sq  is  set  to  a  real  constant,  and  sq  is  called  the  point  spread  function. 
Therefore,  a  conventional  2D  SAR  imaging  operator  is  [76] 

I  POO  PIT 

b{}=^j  oJ  {'K2‘ir|m*/.  (io) 

where  r  is  a  vector  from  the  radar  phase  center  to  the  spatial  coordinates  in  a  chosen 
imaging  plane  at  any  azimuth  angle.  Because  most  of  the  energy  in  a  SAR  scene 
is  well-modeled  by  ideal  point  scatterers,  imaging  operators  which  assume  an  ideal 
point  scatterer,  such  as  Equation  (10),  are  common  and  used  exclusively  throughout 
this  dissertation. 

Note  that  for  2D  imaging,  a  scatterer  does  not  need  to  physically  lie  in  the  imag¬ 
ing  plane  in  order  to  integrate  coherently.  Coherent  integration  occurs  when  |r?|  of 
the  scatterer  equals  |r|  of  the  imaging  operator.  Thus,  scatterers  with  height,  zq,  per- 
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pendicular  to  the  2D  imaging  plane  are  modeled  as  appearing  at  off-set  coordinates 
determined  by  a  projection  along  the  spherical  wavefront  into  the  imaging  plane. 
This  is  effect  is  referred  to  as  SAR  layover  [76].  As  a  result,  |r9|  can  be  modeled  as 
a  function  of  xq,  yq ,  and  9,  where  (. xq,yq )  are  the  offset  coordinates  due  to  layover 
effect,  if  any.  As  such,  the  shorthand  notation 

S,(f,g)e-^  =  S,(f,g-,x„y,)  (11) 

will  be  used  throughout  this  dissertation. 

2. 1.3.1  Observations. 

It  has  been  observed  that  certain  objects  in  a  scene  may  exhibit  non-ideal  scat¬ 
tering  behavior  with  amplitude  functions  that  are  anisotropic,  dispersive,  or  both.  In 
addition,  moving  objects  have  coordinates,  (xq,  yq),  that  vary  with  time.  As  a  result  of 
their  non-ideal  behavior,  such  objects  may  appear  unfocused  or  displaced  in  the  SAR 
image  reconstructed  from  Equation  (10).  The  insertion  of  additional  filters  into  the 
imaging  operator  can  cause  an  anisotropic,  dispersive,  or  moving  object  to  simulate 
ideal  scattering  behavior  and  become  better  focused  in  a  SAR  image,  and  adaptive 
filters,  such  as  those  used  in  multiple  signal  classification  (MUSIC),  can  enhance  the 
detection  and  of  a  pre- determined  type  of  scatterer  [39,  66,  70,  94,  126,  130,  145]. 
However,  such  filters  have  the  negative  consequence  of  defocusing  other  scatterers  of 
interest,  especially  if  the  filters  are  non-adaptive. 

The  imaging  algorithms  with  a  computational  complexity  of  0(N2  log  N)  use  an 
inverse  fast  Fourier  transform  (IFFT)  in  two  dimensions,  where  the  computational 
complexity  of  a  single  IFFT  is  0(N  log  N).  Unfortunately,  an  FFT  for  polar  coordi¬ 
nates  is  not  known.  Therefore,  these  algorithms  transform  the  SAR  phase  history  and 
imaging  algorithm  into  a  rectangular  format  before  using  the  IFFT.  This  transfor- 
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mation  and  its  associated  approximations  produce  many  of  the  limitations  associated 
with  the  faster  imaging  algorithms.  Therefore,  the  characterization  and  mitigation  of 
these  limitations  comprise  much  of  the  literature  on  SAR  imaging. 

This  dissertation  assumes  the  data  and  imaging  operator  are  in  polar  coordinates, 
as  in  Equations  (7)  and  (10),  throughout  its  development.  However,  the  theory  and 
algorithm  are  not  limited  to  polar  coordinates,  and  considerations  for  rectangular 
formatted  data  and  imaging  operators  are  discussed,  where  appropriate. 

2.2  Physical  Model  of  Canonical  Scatterers 

Unfortunately,  the  number  of  ideal  point  scatterers  required  to  accurately  model 
or  simulate  a  SAR  phase  history  is  typically  quite  large,  particularly  for  wide  band 
or  wide  apertures  data.  A  parsimonious  sum  is  possible  when  the  scattering  cen¬ 
ters  are  modeled  as  canonical  scatterers  [62],  Geometric  objects  in  the  scene,  called 
canonical  scatterers,  have  a  predictable  response  to  changes  in  frequency,  azimuth, 
and  polarimetry  [83].  Examples  of  canonical  scatterers  include  trihedrals,  dihedrals, 
plates,  cylinders,  and  spheres.  Canonical  scatterers  are  of  interest  because  they  are 
commonly  associated  with  man-made  objects  [14,  106,  108,  149].  For  example,  the 
Sandia  Laboratories  implementation  of  cylinders  (SLICY)  is  comprised  of  canonical 
scatterers  as  shown  in  Figure  7.  Here,  the  cylinders  are  considered  a  special  case  of 
the  general  cone  shape. 

It  is  desirable  to  detect  and  classify  canonical  scatterers  in  SAR  imagery  using  well- 
known,  physical  models  of  electromagnetic  scattering  [29,  122],  It  has  been  shown  that 
the  amplitude  functions  for  canonical  scatterers  are  parameterized  by  physical  models 
based  on  geometric  optics  (GO)  and  the  Geometric  Theory  of  Diffraction  (GTD)  [78]. 
A  restricted  set  of  possible  geometrical  shapes,  combined  with  the  high-frequency, 
far-held  assumptions  in  GO/GTD,  produce  a  model  with  only  a  few  parameters. 
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Figure  7.  The  SLICY  can  be  modeled  as  a  collection  of  canonical  scatterers  [149].  Note 
that  a  cylinder  can  be  represented  as  a  special  case  of  a  cone. 


The  model  can  be  used  to  simulate  a  SAR  phase  history,  or  the  parameters  can  be 
estimated  from  SAR  data  to  detect  or  classify  canonical  scatterers  in  the  scene. 

In  order  to  develop  a  canonical  scatterer  classification  algorithm  in  the  spatial 
domain,  subimages  are  modeled  as  a  sum  of  subimages  due  to  Q  canonical  scatterers. 
This  idea  is  developed  in  detail  later  in  Chapter  III.  Meanwhile,  the  following  subsec¬ 
tions  present  the  background  needed  to  understand  how  the  amplitude  functions  of 
canonical  scatterers  vary  with  changes  in  frequency  and  azimuth.  They  also  present 
how  the  intensity  of  canonical  scatterers  responds  to  changes  in  polarimetry. 
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cr,rf  =  f-\2jtb'lc  o„,  =  final1 /c  alph  » nata2 

Figure  8.  Canonical  scatterer  radar  cross  sections  have  a  frequency  dependency  that 
depends  upon  the  local  curvature  of  the  scatterer  [83]. 

2.2.1  Amplitude  Response  in  Frequency. 

For  a  co-polarized  channel,  the  amplitude  response  of  the  electric  held  backscat- 
tered  from  a  canonical  scatterer  has  a  frequency  response  predicted  by  GO/GTD  as 
[62,  78] 

Sf(f;A,a)=A(jfr/2,  (12) 

where  A  is  a  complex-valued  amplitude  related  to  the  physical  size  of  the  canonical 
scatterer,  /  is  the  frequency  of  the  incident  electromagnetic  held,  and  a  is  an  integer 
value  depending  upon  the  local  curvature  of  the  canonical  scatterer’s  shape.  Note 
that  this  model  assumes  the  canonical  scatterers  are  perfect  electrical  conductors. 

Figure  8  illustrates  how  frequency  dependency  of  the  radar  cross  section  (RCS) 
for  canonical  scatterers  depends  upon  the  scatterer’s  local  curvature.  In  this  case,  the 
RCS  is  the  magnitude  squared  of  the  amplitude  function  given  by  a  =  \Sf\2  oc  /“. 
The  trihedral,  having  no  curvature,  features  a  quadratic  response  (/2);  the  cylinder, 
having  curvature  in  one  dimension,  features  a  linear  response  (f1)’,  and  the  sphere, 
having  curvature  in  two-dimensions,  features  a  hat  response  (/°).  The  values  of  a 
for  common  shapes  are  well  known  and  listed  in  Table  2.  Note  that  there  is  an 
ambiguity  when  discriminating  between  canonical  scatterers  by  a  only.  For  instance, 
plates,  trihedrals,  and  dihedrals  all  have  the  frequency  parameter,  a  =  2. 
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Table  2.  Traditional  Frequency  Response  Parameter  for  Ideal  Canonical  Scattering 
Geometries  [83]. 


Scattering  geometry 

OL 

plate,  trihedral,  dihedral 

2 

cylinder /cone,  top  hat 

1 

sphere,  straight  edge/wire 

0 

Note  that  diffraction  from  curved  edges  and  tips  of  objects  produce  scattering  with 
inverse  frequency  response  of  a  =  —1  and  a  =  —2,  respectively  [62,  84,  121].  However, 
these  have  such  a  low  RCS  as  to  not  be  prevalent  in  SAR  imagery.  As  a  result,  curved 
edge  and  tip  scattering  mechanisms  are  ignored  throughout  this  dissertation,  except 
in  Section  5.4.2. 

2.2. 1.1  Observations. 

The  GO/GTD-based  scattering  models  are  limited  by  a  high-frequency  approxi¬ 
mation.  However,  these  models  are  preferred  because  backscattering  which  occurs  at 
lower  frequencies  has  less  directivity,  which  lowers  the  received  energy  of  the  desired 
signal.  For  wavelengths  greater  than  the  object  extent,  the  frequency  response  is  gov¬ 
erned  by  Rayleigh  scattering  [83].  In  this  case,  there  is  little  variation  of  the  incident 
held  across  the  object,  and  the  incident  held  can  be  modeled  as  being  quasi-static 
[83].  Under  these  conditions,  relative  intensities  of  scatterers  can  be  determined,  but 
the  amplitude  response  is  independent  of  object  shape  [83].  Because  shape  cannot  be 
determined,  no  models  exist  for  describing  canonical  scatterers  when  the  wavelength 
of  the  incident  held  is  greater  than  the  extent  of  the  canonical  scatterer. 

Alternately,  when  wavelengths  approach  the  order  of  the  object  size,  it  has  been 
shown  that  an  object’s  size  is  related  to  its  late-time  resonance  response  [27,  79].  In 
this  case,  physical  mechanisms  cause  EM  energy  to  stay  attached  to  the  surface  of  an 
object  in  what  are  called  surface  waves  [83].  Surface  wave  scattering  produces  a  strong 
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early-time  signal  that  is  independent  of  object  size  [83].  However,  it  has  been  shown 
that  the  much  weaker  late-time  signal  will  resonate  in  accordance  with  object  size 
[79].  Models  describing  this  late-time  resonance  response  have  been  successfully  used 
to  classify  objects,  but  only  in  ideal  scenarios  where  noise,  clutter,  and  interference 
are  constrained. 

In  general,  the  performance  of  late-time  resonance  approaches  to  scatterer  clas¬ 
sification  is  limited  because  of  the  low  energy  of  the  late-time  response.  Thus,  in 
practice,  accurate  feature  extraction  requires  a  prohibitively  high  signal-to-noise  ra¬ 
tio.  Matched  filter  techniques,  such  as  those  found  in  singularity  expansion  method 
[9] ,  have  been  shown  to  reduce  noise  sensitivity,  but  at  the  expense  of  requiring  a  priori 
knowledge  of  the  target.  Despite  these  limitations,  interest  in  late-time  resonance  re¬ 
sponse  models  continues,  as  evidenced  by  recent  publications  [61,  105].  Nevertheless, 
limiting  factors  restrict  the  usefulness  of  resonance  response  models  in  distinguishing 
canonical  scatterers. 

2.2.2  Amplitude  Response  in  Azimuth. 

At  a  single  azimuth  angle  or  over  an  extremely  narrow  aperture,  the  amplitude 
responses  for  all  canonical  scatterers  are  well  modeled  by  Equation  (12)  [121],  How¬ 
ever,  for  typical  SAR  apertures,  the  amplitude  response  of  the  baekscattered  held 
for  distributed  scatterers  has  an  azimuth  dependency  dominated  by  a  sinc-like  pat¬ 
tern.  This  response  is  in  relation  to  the  slant  plane  containing  the  synthetic  aperture. 
Common  distributed  scatterers  include  hat  plates  at  broadside  aspect,  dihedrals  with 
fold-lines  parallel  to  the  slant  plane,  cylinders  with  axes  of  rotation  parallel  to  the 
slant  plane,  and  edges  or  wires  lying  parallel  to  the  slant  plane.  Of  these,  the  di¬ 
hedral,  in  particular,  is  often  present  in  SAR  imagery  of  man-made  structures.  For 
example,  the  side  of  a  building  and  the  ground  form  a  dihedral  with  a  fold-line  often 
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synthetic  aperture 


Figure  9.  Example  model  parameters  for  a  cylinder  with  axis  of  rotation  projected  into 
the  slant  plane. 


lying  parallel  to  the  slant  plane.  The  amplitude  response  as  a  function  of  azimuth  is 
well-modeled  by  [6,  62] 


Se(f,  0]  L,  d0)  =  sine  [~cfL  sin(d  -  d0)]  ,  (13) 

where  the  sine  function  is  defined  as  sinc(i)  =  ,  L  is  the  effective  length  of  the 

scatterer  as  projected  onto  the  slant  plane,  and  do  is  the  orientation  angle  normal 
to  this  projection  and  referenced  to  the  center  angle  of  the  aperture,  dc  [62].  The 
dependence  upon  dc  is  suppressed  in  the  notation  because  the  simplification  dc  =  0 
can  often  be  made  without  loss  of  generality.  An  example  of  the  effective  length  and 
orientation  angle  for  a  cylinder  is  illustrated  in  Figure  9. 
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It  has  been  shown  that  most  scatterers  have  an  angular  persistence  of  less  than 
twenty  degrees  [43].  Furthermore,  for  flat  plates,  there  is  an  extra  cos(0)  factor  in 
the  azimuth  response  of  Equation  (13)  [6].  For  these  reasons,  it  is  prudent  to  place 
a  restriction  of  0  <  20°  on  the  aperture  width  for  which  Equation  (13)  is  valid. 
However,  in  practice  any  aperture  greater  than  20°  will  likely  be  decomposed  into 
narrower  subapertures  before  imaging;  so  this  restriction  is  deemphasized  throughout 
the  dissertation. 

2. 2. 2.1  Observations. 

For  very  narrow  angle  imaging,  scatterer  anisotropy  is  generally  negligible  because 
large  variations  in  target  aspect  are  not  expected  [125].  However,  anisotropy  should 
be  accounted  for  any  time  large  variations  in  aspect  are  encountered,  such  as  for  wide- 
angle  SAR  applications  [134],  These  include  strip-map  data  collected  from  air  and 
spaceborne  radar  platforms  operating  at  P-  or  L-band  [51],  data  collected  from  RCS 
measurement  facilities  where  targets  are  placed  on  turn-tables,  and  circular  SAR  data 
collected  from  airborne  radar  platforms  [110].  Much  of  the  early  research  in  target 
anisotropy  was  motivated  by  radar  imaging  at  low-frequencies  for  foliage  and  ground 
penetration.  The  use  of  low- frequencies  drove  this  early  research  to  examine  the 
late-time  resonance  response  of  targets  [8,  85]. 

Before  the  mid-1990s,  such  wide-angle  SAR  data  collections  were  not  common. 
However,  recent  technological  advances  have  enabled  the  collection  of  coherent  data 
over  very  wide  apertures  [110].  These  developments  have  spurred  interest  over  the 
last  decade  in  studying  methods  and  models  that  predict  and  leverage  the  anisotropic 
behavior  of  scatterers.  A  review  of  the  current  literature  reveals  that  directional  filters 
are  the  most  common  approach  to  leveraging  anisotropic  behavior  for  target  detection 
and  discrimination  [1,  47,  130,  136].  Directional  filters  denote  subaperture  techniques 
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where  the  subaperture  filter  is  often  chosen  according  to  expected  target  anisotropy. 
When  ultra-wide  band  data  is  available,  the  directional  filters  are  often  extended  to  2D 
filtering  to  account  for  both  anisotropic  and  dispersive  analysis.  In  addition  to  filters, 
the  parametric  model  in  Equation  (13)  [31,  50,  74]  and  a  sparse  dictionary  method 
[144]  have  also  been  proposed  to  characterize  and  identify  anisotropic  scattering. 

2.2.3  Intensity  Response  in  Polarization. 

This  subsection  presents  the  theoretical  background  needed  to  understand  how  the 
intensities  of  canonical  scatterers  are  polarization  dependent.  Specifically,  it  presents 
Krogager  decomposition  of  the  scattering  matrix  as  a  way  to  differentiate  between 
odd-bounce  and  even-bounce  scattering  mechanisms.  The  proportions  of  odd-bounce 
and  even-bounce  scattering  energy  are  very  useful  in  determining  the  geometric  shape 
of  canonical  scatterers. 


2.2.3. 1  The  Scattering  Matrix. 


The  scattering  matrix  can  be  used  to  specify  polarimetric  properties  of  electromag¬ 
netic  scattering.  That  is,  the  polarization  pairs  of  the  scattered  held  are  determined 
via  matrix  multiplication  of  the  scattering  matrix  with  the  polarization  pairs  of  the 
incident  held  [112] 
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where  the  elements  of  the  scattering  matrix  are  complex,  e  C,  and  r  is  the  distance 

between  the  receive  antenna  and  the  reference  plane  at  the  scatterer.  The  polarization 
pairs  can  be  expressed  in  either  linear  or  circular  polarizations.  This  choice  is  often 
dictated  by  the  antenna  design  of  the  radar  system.  Standard  coordinate  system 
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conventions  for  linear  polarization  pairs  use  x  =  H  for  horizontal  polarization  and  y  = 
V  for  vertical  polarization  [112].  While,  similar  conventions  for  circular  polarization 
pairs  use  x  =  L  for  left  circular  and  y  =  R  for  right  circular  [112]. 


The  definition  of  the  scattering  matrix  is  not  unique  because  the  existence  of 
many  equivalent  coordinate  systems  introduces  an  ambiguity.  However,  two  particular 
coordinate  systems  and  resulting  scattering  matrix  definitions  are  commonly  used. 
The  Jones  matrix  is  popular  because  it  presents  a  right-handed  system  with  regard 
to  the  conventional  definition  of  wave  propagation  [112].  In  contrast,  the  Sinclair 
matrix,  S,  presents  a  left-handed  system.  A  left-handed  system  would  not  normally 
be  desirable,  but  the  Sinclair  matrix  has  the  distinct  advantage  that  Sxy  =  Syx  for 
the  case  of  monostatic  backscattering.  The  dissertation  is  limited  to  the  case  of 
monostatic  backscattering  only;  therefore,  the  Sinclair  scattering  matrix  given  by 

(15) 

will  be  used  exclusively  throughout  the  remaining  discussion.  In  this  case,  the  trans¬ 
formation  from  the  linear  polarization  basis  to  the  circular  polarization  basis  is  given 
by  [88] 

Srr  =  JSrv  +  \{Shh  —  SVv) 

Sll  =  jSnv  ~  \{Shh  —  Svv)  (16) 

Srl  =  1(Shh  +  Svv), 
where  Slr  =  Srl  for  monostatic  radar  systems. 

2. 2. 3. 2  Krogager  Decomposition  of  the  Sinclair  Scattering  Matrix. 

Several  decompositions  of  the  Sinclair  scattering  matrix  have  been  proposed  for 
canonical  scatterer  analysis.  The  most  common  being  Pauli  [33],  Krogager  [87],  and 
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Figure  10.  Reflection  behavior  for  linearly  polarized  electric  fields.  The  orientation  of 
both  the  vertical  and  horizontal  components  of  the  electric  field  are  unperturbed  after 
an  odd  number  of  bounces  (top  and  bottom).  However,  the  vertical  component  of  the 
electric  field  becomes  inverted  following  a  double  bounce  (middle)  [6]. 


Cameron  [20].  Figure  10  uses  linear  polarization  pairs  to  illustrate  the  motivation 
behind  such  canonical  scatterer  analysis.  The  orientations  of  both  the  vertical  and 
horizontal  components  of  the  electric  held  are  unperturbed  after  a  single  bounce  as 
shown  in  the  top  diagram.  In  contrast,  the  vertical  component  of  the  electric  held, 
Ey,  becomes  inverted  after  a  double  bounce  as  shown  in  the  bottom  diagram.  After 
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a  triple  bounce,  such  as  may  be  experienced  in  the  well  of  a  trihedral,  a  double 
inversion  will  cause  the  scattered  held  to  return  to  its  original  orientation  as  shown  in 
the  bottom  diagram.  As  such,  the  vertical  component  of  the  electric  held  is  inverted 
for  all  even-bounce  geometries  and  is  unchanged  for  all  odd-bonnce  geometries.  The 
linear  polarization  Sinclair  matrix  for  ideal  odd-bounce  geometries  is  thus  given  by 
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1  0 
0  1 


(17) 


While  the  linear  polarization  Sinclair  matrix  for  ideal  even-bounce  geometries  is  given 
by 


Se 


1  0 
0  -1 


(18) 


In  fact,  these  matrices  comprise  two  terms  of  the  matrix  decomposition  based  on 
Pauli  spin  matrices  (referred  to  as  Pauli  decomposition)  [92] 
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where  a,  b ,  c,  and  d  are  complex  valued  and  given  by  [92] 


Shh  +  Sw  ,  Shh  —  Syv  Shv  +  Svh  ,  .Shv  —  Svh 

a  = -  b  = -  c  = -  d  =  / - .  (20) 
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Note  that  for  monostatic  backscatter,  d  =  0. 

In  order  to  better  understand  Pauli  decomposition,  it  is  benehcial  to  further  ex¬ 
amine  the  double-bounce  mechanism  depicted  in  the  middle  diagram  of  Figure  10. 
Note  that  the  horizontal  component  of  the  incident  electric  held,  ElH,  is  parallel  to 
the  dihedral  fold  line.  If  the  situation  were  altered  by  rotating  the  dihedral  fold  line 


30 


by  90°,  then  the  vertical  component  of  the  incident  electric  held,  Ey,  would  become 
parallel  to  the  fold  line.  Conversely  and  according  to  this  new  orientation,  Eh  would 
be  inverted  while  Ey  would  remain  unchanged  overall.  Such  a  change  in  orientation 
is  referred  to  as  a  roll.  In  Equation  (19),  a  90°  roll  would  simply  invert  the  sign  of  b, 
leaving  the  magnitude  unchanged.  However,  in  general,  the  magnitude  of  b  depends 
upon  roll  angle. 

The  third  term  of  Equation  (19)  can  be  interpreted  to  represent  backscatter  from 
a  dihedral  with  a  fold  line  rotated  45-degrees,  as  compared  to  that  shown  in  the 
middle  diagram  of  Figure  10  [92],  Therefore,  dihedrals  of  arbitrary  roll  contribute  to 
both  the  second  and  third  terms  of  the  Pauli  decomposition.  In  response,  Krogager 
proposed  a  roll-invariant  decomposition  of  the  Sinclair  scattering  matrix  [86,  87].  The 
Krogager  decomposition  expressed  in  circular  polarization  bases  can  be  modeled  as 
[91] 

Stir  Srl 
Slr  Sll 

where  the  phase  difference  between  HH  and  VV  polarizations  is  assumed  to  be  zero 
and  (p  is  called  the  helix  phase  angle. 

The  coefficients  K01  Ke,  and  Kh  can  be  extracted  from  polarization  diverse  SAR 
images  on  a  pixel-by-pixel  basis  using  [90] 


10  10  0  1 

=  K0  +Ke  +Khe ,  (21) 
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K0  =  |Srl|, 

Ke  =  minds'll,  \SRR\),  (22) 

Kh  =  abs(|S'_R,i?|  —  |Sll|)|, 

where  the  coefficients  represent  the  strength  of  odd-bounce,  even-bounce,  and  helical 
scattering,  respectively.  This  interpretation  is  valid  for  all  canonical  scatterers  [92], 
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Table  3.  Krogager  and  Frequency  Response  Parameters  for  Canonical  Point  Scattering 
Geometries. 


Scattering  geometry 

at. 

K-q 

flat  plate  or  trihedral 

2 

1 

0 

dihedral 

2 

0 

l 

cylinder/cone  specular 

1 

l 

0 

top  hat 

1 

0 

1 

sphere 

0 

l 

0 

straight  edge/wire  specular 

0 

0.5 

0.5 

The  relative  strengths  of  odd-bounce  and  even-bounce  scattering  can  be  deter¬ 
mined  by 
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Ke  | 

|2  + 

\Kh\ 2 

y/\K.\ 

P  + 

|  Ke 

|2  + 

\Kh\* 

(23) 

(24) 


where  k0  and  Ke  are  real-valued,  normalized  scattering  intensities  measured  at  a 
given  image  location  or  pixel.  In  this  case,  the  amount  of  helical  scattering  can  be 
interpreted  as  relating  to  the  purity  of  the  even-bounce  scattering  [88]. 

By  simply  combining  these  Krogager  parameters  with  the  frequency  parameter 
from  Table  2,  some  of  the  ambiguity  between  scatterers  can  be  resolved  in  three  di¬ 
mensions.  A  listing  of  ideal  canonical  point  scatterers  based  on  these  three  parameters 
is  given  in  Table  3. 


2. 2. 3. 3  Observations. 

The  goal  of  scattering  matrix  decomposition  is  to  produce  a  set  of  basis  matrices 
which  can  give  insight  into  the  type  of  scattering  present  in  a  radar  signal  [92] .  Such 
insight  has  been  shown  to  aid  target  detection  and  image  segmentation  as  described 
in  the  survey  books  [92,  112]  and  papers  [33,  140].  These  surveys  describe  how 
scattering  matrix  decomposition  is  fundamentally  driven  by  the  time-varying  nature 
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of  the  target’s  polarimetric  signal.  If  the  signal  can  be  approximated  as  time-invariant, 
then  coherent  scattering  matrix  decompositions  are  appropriate.  However,  if  the 
signal  is  time- variant,  then  power-type  matrix  decompositions  are  required  [112]. 

When  applicable,  coherent  scattering  matrix  decompositions  are  preferred  because 
they  can  better  describe  physical  scattering  mechanisms  [112].  Furthermore,  coherent 
scattering  matrices  for  backscattering  have  only  five  independent  variables  as  com¬ 
pared  to  ten  for  the  power-type  matrices  [112].  Fortunately,  the  scatterer  classification 
algorithm  developed  in  this  dissertation  only  extracts  polarimetric  features  from  pix¬ 
els  where  the  response  from  a  single  canonical  scatterer  dominates.  In  this  case, 
the  polarization  response  can  be  approximated  as  time-invariant,  and  the  coherent 
scattering  matrix  decomposition  is  applicable  [92], 

There  are  three  well-established  coherent  scattering  matrix  decompositions  that 
have  been  proposed  in  the  literature.  These  are  Pauli  [33],  Krogager  [87],  and 
Cameron  [20].  The  Pauli  decomposition  is  adopted  from  the  Pauli  spin  matrices 
originating  in  physics  and  optics.  It  results  in  orthogonal  basis  matrices  which  for 
radar  backscattering  can  be  interpreted  as  a  sphere,  a  dihedral  with  a  zero-degree  roll, 
and  a  dihedral  with  a  45-degree  roll.  The  Krogager  decomposition  is  presented  as  a 
roll  invariant  alternative  to  the  Pauli  decomposition.  It  has  the  advantage  of  separat¬ 
ing  odd-bounce  from  even-bounce  backscattering.  The  basis  matrices  for  Krogager 
decomposition  can  be  interpreted  as  a  sphere  (odd-bounce),  a  dihedral  (even-bounce), 
and  a  helix.  A  disadvantage  of  Krogager  decomposition  is  that  the  basis  matrices  are 
not  orthogonal.  This  results  in  backscatter  from  double  dihedrals  to  appear  as  helical. 
The  Cameron  decomposition  attaches  a  great  deal  of  importance  to  a  class  of  targets 
which  Huynen  termed  symmetric  [71].  A  symmetric  target  has  an  axis  of  symmetry 
in  the  plane  orthogonal  to  the  radar  line  of  site.  For  these  targets,  decomposition 
follows  a  decision  tree  where  very  detailed  target  information  can  be  gleaned. 
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Table  4.  Comparison  of  Coherent  Scattering  Matrix  Decompositions. 


Decomp. 

Recent  Publications 

Advantages 

Limitations 

Pauli 

2010  Inci  [72] 

2010  Brigui  [16] 

orthogonal  basis 

not  roll  invariant 

Krogager 

2008  Ainsworth  [2] 

2010  Zou  [154] 

roll  invariant 

not  an  orthogonal  basis 

Cameron 

2009  Cameron  [21] 

2009  Martorella  [101] 

details  symmetric  targets 

wide  apertures  degrade  symmetry 
high-order  feature  space 

Touzi 

2008  Brisco  [17] 

roll  invariant 
orthogonal  basis 

focused  toward  non-coherent  decomposition 
new  and  relatively  unproven 

It  is  important  to  note  that  Touzi  has  also  made  some  significant  contributions 
to  the  study  of  coherent  scattering  matrix  decompositions.  He  recently  presented 
a  coherent  decomposition  that  is  both  orthogonal  and  roll  invariant  [139],  combin¬ 
ing  the  strengths  of  the  Pauli  and  Krogager  decompositions,  respectively.  His  intent 
was  to  produce  a  decomposition  that  was  useful  for  both  coherent  scattering  matrix 
and  power-type  matrix  decomposition.  Cameron  recently  raised  the  question  of  the 
uniqueness  of  his  approach  [21],  and  at  this  time,  it  is  unclear  if  the  Touzi  decomposi¬ 
tion  will  rival  the  popularity  of  the  three  well-established  coherent  scattering  matrix 
decompositions  previously  discussed. 

A  summary  of  the  advantages  and  limitations  of  each  method  is  provided  in  Ta¬ 
ble  4.  Because  the  primary  contribution  of  this  dissertation  lies  in  the  development 
of  a  new  method  for  characterizing  the  anisotropic  and  dispersive  characteristics  of 
scatterers,  the  choice  of  scattering  matrix  decomposition  method  is  secondary.  The 
Krogager  decomposition  is  chosen  because  of  its  simplicity  and  roll-invariance.  Al¬ 
though  the  Cameron  decomposition  would  potentially  provide  more  information,  this 
comes  at  the  expense  of  increasing  the  dimensions  of  the  feature  space  and  at  the 
restriction  of  narrow  apertures  to  preserve  symmetry. 
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2.3  Feature  Extraction  and  Classification 


This  section  provides  an  introduction  to  some  basic  concepts  in  feature  extraction 
and  classification,  as  well  as  an  overview  of  the  feature  extraction  and  least  squares 
scatterer  classification  method  featured  in  this  dissertation.  The  subsections  present 
some  specific  examples  of  feature  selection,  feature  extraction,  and  least  squares  clas¬ 
sification  for  canonical  scatterers.  These  focus  on  use  of  the  parameters  for  frequency, 
azimuth,  and  polarization  dependency  discussed  previously  in  this  chapter. 

2.3.1  Feature  Selection. 

Feature  selection  is  an  engineering  art,  where  the  goal  is  to  replace  unwieldy 
high-order  representations  of  objects  with  elegant  lower-order  approximations  [108]. 
Within  the  limits  of  the  approximation,  this  allows  automated  object  detection  and 
classification  to  be  tractable  on  modern  computers  using  signal  detection  and  estima¬ 
tion  theory.  The  features  are  derived  from  model-based  parameters,  statistical-based 
parameters,  or  both.  An  example  model-based  parameter  is  the  order  of  the  amplitude 
response  with  frequency,  parameterized  by  a  in  Equation  (12)  of  Section  2.2.1,  which 
is  shown  in  Table  3  to  be  useful  in  discriminating  between  some  types  of  scatterers. 

An  example  statistical-based  parameter  is  the  expected  value  of  the  pixel  intensity 
associated  with  a  given  object.  For  instance,  at  X-band  (/  =  9  GHz)  and  referring  to 
Figure  8  of  Section  2.2.1,  the  radius  of  a  sphere  must  be  360  meters  in  order  to  produce 
the  same  backscattered  energy  of  a  1  meter  long  trihedral.  Because  such  a  large 
sphere  is  not  expected  in  the  scene,  certain  high- intensity  pixels  are  not  expected  to 
be  associated  with  a  sphere.  Alternately,  statistical  parameters  are  often  derived  from 
statistical  pattern  recognition  techniques  using  training  data.  An  example  of  training 
data  is  phase  histories  from  scenes  containing  known  objects  at  known  locations. 
Because  training  data  is  scenario  specific  and  this  dissertation  seeks  a  general  theory 
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for  scatterer  classification,  statistical  parameters  are  excluded  from  the  remaining 
discussion. 

2.3.2  Feature  Extraction. 

Model-based  feature  extraction  is  accomplished  by  a  series  of  linear  and  non-linear 
estimation  techniques,  depending  upon  the  types  of  features  being  extracted.  For 
instance,  model  parameters  are  easily  extracted  from  SAR  data  by  fitting  measured 
data  points  and  curves  to  responses  predicted  by  the  models.  An  example  is  image 
peak  detection,  which  can  be  used  to  estimate  the  location  (xq,  yq)  of  the  qth  scattering 
center.  Another  example  is  least  squares  fitting  of  the  main  lobe  of  the  sine  function 
in  azimuth  to  a  quadratic  approximation  in  the  spectral  domain  parameterized  by 
L  [3].  A  final  example  is  pixel  summation  from  images  of  different  polarizations  to 
obtain  the  Krogager  coefficients. 

After  the  parameters  are  extracted,  they  are  recorded  in  a  feature  vector.  For 
canonical  scatterers,  all  of  the  model  parameters  previously  introduced  in  this  chapter 
can  be  arranged  in  the  feature  vector 

w  =  [; x,y ,  \A\,ZA,a,L,90,Ko,Ke].  (25) 

2.3.3  Least  Squares  Scatterer  Classification. 

Scatterer  classification  is  accomplished  by  a  likelihood  ratio  test,  where  the  mea¬ 
sured  feature  vector  is  compared  to  a  set  of  ideal  feature  vectors,  each  representing  a 
particular  scatterer  type.  The  simplest  form  of  likelihood  ratio  test  is  the  least  squares 
classifier,  where  the  vectors  are  represented  in  a  Euclidean  space  and  are  compared 
using  the  Euclidean  norm.  The  simplicity  of  the  least  squares  classifier  makes  it  the 
clearest  method  by  which  to  test  the  hypotheses  of  this  dissertation.  Once  the  new 
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Figure  11.  An  example  2D  feature  space  for  a  least  squares  classifier  and  four  classes. 
The  extracted  feature  vector,  vq,  is  superimposed  on  a  feature  space  divided  into  four 
regions  by  classification  basis  represented  by  ideal  feature  vectors  weven,  w edgei  w odd, 
and  w heUx. 


theory  for  scatterer  classification  via  domain  decomposition  is  established,  follow-on 
work  with  more  sophisticated  classifiers  can  be  accomplished. 

An  illustration  of  least  squares  classification  with  a  two-dimensional  feature  vec¬ 
tor,  w  =  [k,0,  Ke],  is  given  in  Figure  11.  There  are  four  ideal  feature  vectors,  each 
corresponding  to  a  different  class  of  scatterer.  The  Euclidean  norm  divides  the  fea¬ 
ture  space  into  four  regions,  one  for  each  class.  An  example  extracted  feature  vector 
for  the  qth  scatterer  is  illustrated  by  vq,  which  is  nearest  weven,  as  measured  by  the 
Euclidean  norm.  In  this  case,  the  gth  scatterer  is  most  likely  a  canonical  scatterer 
producing  an  even-bounce,  such  as  a  dihedral. 

2.3.3. 1  Observations. 

The  entire  gamut  of  pattern  recognition  techniques  have  been  adapted  for  char¬ 
acterizing  and  identifying  objects  in  SAR  imagery  [63,  142],  While,  some  of  these 
techniques  have  proven  useful  in  different  applications,  many  approaches  are  limited 
by  the  dilemma  of  combinatorial  complexity,  especially  as  the  number  of  features 
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increases  [119].  In  contrast,  the  use  of  simple  physical  models  with  only  a  few  param¬ 
eters  has  the  advantage  of  producing  efficient  and  optimal  solutions,  by  the  principle 
of  Occam’s  razor  [11], 

Furthermore,  scatterer  classification  results  can  be  used  as  a  first  step  in  a  tree 
algorithm  or  as  an  initialization  mechanism  for  higher-order  classification  algorithms. 
For  example,  Reference  [149]  describes  how  information  gleaned  from  a  canonical  scat¬ 
terer  classifier  enables  a  secondary  target  classification  algorithm  to  produce  better 
separation  between  classes  for  subsequent  hypothesis  testing. 

A  least  squares  classifier  is  optimal  when  the  desired  signal  is  corrupted  only 
by  additive  white  Gaussian  noise  (AWGN).  The  thermal  noise  of  a  SAR  system 
produces  AWGN  in  the  phase  history  [133].  Likewise,  the  image  contains  AWGN  with 
a  scaled  variance  clue  to  coherent  processing  of  the  imaging  operator  [133].  However, 
because  the  frequency  and  polarimetric  parameters  are  taken  from  the  magnitude  or 
intensity  of  an  image,  the  noise  affecting  these  parameters  is  expected  to  be  colored. 
Furthermore,  the  interference  due  to  neighboring  scatterers  and  clutter  is  expected 
to  be  colored  as  well.  As  a  result,  the  performance  of  the  least  squares  classifier  used 
in  this  dissertation  is  expected  to  be  suboptimal. 

Fortunately,  there  are  some  scenarios  where  the  noise,  clutter,  and  interference 
can  be  minimized,  particularly  for  simple  targets  in  free-space.  Examples  include 
stealthy  aircraft  in  flight,  streamlined  spacecraft,  and  simple  objects  inside  anechoic 
chambers.  For  other  scenarios,  such  as  imaging  of  the  earth’s  surface  or  imaging  of 
complex  targets  comprised  of  many  scatterers  in  close  proximity,  the  performance  of 
the  least  squares  classifier  is  expected  to  be  suboptimal.  Nonetheless,  a  least  squares 
classifier  is  sufficient  to  illustrate  the  usefulness  of  the  phase  history  decomposition 
method  for  scatterer  classification. 
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III.  Research  Objective 


This  chapter  serves  as  a  foundation  for  the  work  in  this  dissertation.  It  presents 
the  motivation  for  the  research  as  well  as  some  of  the  expected  benefits  and  limitations 
of  the  new  approach  to  imaging  and  scatterer  classification  based  on  phase  history 
decomposition. 

3.1  Motivation 

SAR  images  contain  bright  spots  representing  locations  where  strong  backscat- 
terers  are  present  in  the  scene.  Within  these  bright  spots  are  pixels  having  localized 
peak  intensities,  called  image  peaks.  For  example,  Figure  12  is  a  SAR  image  of  a  resi¬ 
dential  scene  taken  from  the  moving  and  stationary  target  acquisition  and  recognition 
(MSTAR)  data  set,  where  bright  spots  and  image  peaks  are  noticeable  throughout 
the  image  [38,  80]. 

Because  the  imaging  operator  assumes  a  scene  consisting  of  ideal  point  scatterers, 
the  dispersive,  anisotropic,  and  polarimetric  (DAP)  characteristics  of  bright  spots  and 
peaks  are  not  revealed  in  the  image.  Rather,  these  characteristics  must  be  determined 
through  additional  analyses.  By  analyzing  the  DAP  characteristics  of  bright  spots 
and  peaks,  it  is  possible  to  measure  the  likelihood  that  a  bright  spot  corresponds  to  a 
specific  type  of  object.  So,  one  may  ask:  By  analyzing  these  characteristics  in  Figure 
12,  is  it  possible  to  determine  whether  the  bright  spots  in  ovals  A  and  B  are  more 
likely  due  to  automobiles  or  construction  equipment  or  is  it  possible  to  determine  if 
the  bright  line  in  oval  C  is  more  likely  due  to  a  fence  or  pipeline? 

In  general,  the  number  of  possible  inquiries  is  unlimited,  and  reasonable  answers 
depend  upon  contextual  clues  which  are  best  deduced  by  human  operators.  Therefore, 
it  is  desirable  to  succinctly  present  the  DAP  characteristics  of  bright  spots  and  peaks 
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Figure  12.  SAR  image  of  a  residential  scene.  Ovals  A,  B,  and  C  represent  bright  spots 
of  interest. 

within  the  context  of  the  image.  Because  SAR  signals  are  processed  in  both  the  spatial 
and  spectral  domains,  it  is  helpful  to  categorize  dispersive  and  anisotropic  analysis 
methods  into  two  categories  —  image  segmentation  and  phase  history  decomposition 
-  according  to  the  domain  in  which  they  decompose  the  signal.  (The  most  common 
polarimetric  analyses  are  performed  in  the  spatial  domain). 

Image  segmentation  methods  decompose  the  SAR  signal  in  the  spatial  domain. 
The  basic  steps  of  the  image  segmentation  methods  are  depicted  in  Figure  13.  First, 
the  segments  are  localized  to  bright  spots  or  known  locations  of  interest,  such  as 
segments  A,  B,  and  C.  Then,  a  time-frequency  transform  is  applied  to  each  segment 
to  produce  a  coarse  resolution  spectrum  for  further  analysis.  Image  segmentation 
methods  are  best  represented  by  two  research  projects:  the  hyperimage  concept  de¬ 
veloped  at  ONERA,  the  French  Aerospace  Lab,  [46,  47,  48]  and  parametric  scatterer 
classification  developed  at  The  Ohio  State  University  [3,  30,  62,  73].  Other  image 
segmentation  methods  in  the  literature  are  generally  a  variation  of  these  two  primary 
methods. 

The  hyperimage  concept  simply  displays  the  2D  resultant  of  a  time-frequency 
transform.  An  example  hyperimage  is  shown  in  Figure  14.  Unfortunately,  a  human 
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Figure  13.  Basic  steps  of  the  image  segmentation  methods. 

analyst  must  manually  select  the  image  segments  and  transforms.  Then  the  analyst 
must  scrutinize  the  resultants,  each  of  which  are  displayed  separately  and  outside  the 
image.  Thus,  the  hyperimage  concept  is  best  suited  for  highly  trained  analysts  having 
ample  time  to  conduct  the  analyses. 

Parametric  scatterer  classification  goes  the  extra  step  to  match  image  segments 
and  their  spectra  to  the  spectral  and  polarimetric  responses  expected  for  a  set  of  ideal 
canonical  scatterers.  The  set  of  canonical  scatterers  is  limited,  so  that  classification 
results  can  be  succinctly  represented  by  a  small  set  of  symbols  and  displayed  as  an 
overlay  within  the  image  context.  Alternately,  the  classification  results  can  be  used 
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Figure  14.  SAR  image  of  a  helicopter  with  example  hyperimage  displaying  the  spectra 
of  seven  image  segments  [46]. 


to  simulate  the  phase  history  as  a  sum  of  canonical  scatterers.  From  this  simulated 
phase  history,  the  bright  spots  in  the  original  image  can  be  reconstructed  as  shown 
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Figure  15.  SAR  image  of  T-72  tank  from  measured  data  (top)  and  simulated  SAR  image 
reconstructed  from  a  canonical  scatterer  model  (bottom).  The  model  parameters  were 
estimated  using  parametric  scatterer  classification  with  human  supervision  [30,  111]. 


in  Figure  15.  Classification  accuracy  is  highly  dependent  upon  the  non-linear  image 
segmentation  process,  and  for  typical  scenes,  human  supervision  is  needed  to  ensure 
quality  image  segmentation  [3,  62,  73].  In  addition,  when  image  segments  contain 
energy  from  more  than  one  canonical  scatterer,  the  selection  of  model  order,  Q, 
also  requires  human  supervision  to  ensure  estimates  of  the  model  parameters  are 
accurate  [3,  62,  73].  Finally,  classification  accuracy  suffers  when  image  segments 
contain  non-canonical  scatterers,  such  as  resonant  cavities.  Therefore,  in  practice, 
parametric  scatterer  classification  demands  human  supervision  at  some  point  in  the 
process  to  ensure  accuracy.  Unfortunately,  existing  image  segmentation  methods 
cannot  succinctly  present  the  DAP  characteristics  of  bright  spots  and  peaks  within 
the  image  context  in  an  operationally  efficient,  automated  way. 
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Figure  16.  Basic  steps  of  the  phase  history  decomposition  methods. 


Phase  history  decomposition  methods  decompose  the  SAR  signal  in  the  spectral 
domain,  as  previously  described  in  Section  2. 1.2. 2.  The  basic  steps  of  phase  history 
decomposition  methods  are  depicted  in  Figure  16.  The  phase  history  is  subdivided  by 
regular  intervals  into  subdomains,  and  coarse  resolution  subimages  are  reconstructed 
from  each  subdomain  for  further  processing.  In  practice,  angular  bandwidth  is  less 


44 


restricted  than  frequency  bandwidth;  therefore,  the  most  common  decompositions 
are  in  azimuth,  as  previously  described  in  2. 1.2. 4.  The  full  aperture  is  subdivided 
into  two  or  more  subapertures  of  equal  width.  The  resulting  subimages  are  diverse  in 
azimuth,  but  with  coarser  resolution  in  cross-range  than  a  full  aperture  image.  These 
subimages  reveal  the  anisotropic  scattering  behavior  of  each  image  pixel,  which  can 
be  analyzed,  or  in  the  case  of  multilook  SAR,  averaged  over  the  full  aperture. 

Similarly,  for  decompositions  in  frequency,  the  bandwidth  is  subdivided  into  two 
or  more  subbands.  Because  the  available  bandwidth  is  often  limited,  the  subbands 
are  usually  equal  to  half  of  the  full  bandwidth.  The  resulting  subimages  are  diverse 
in  frequency,  but  with  coarser  resolution  in  range  than  a  fullband  image.  These 
subimages  reveal  the  dispersive  scattering  behavior  of  each  image  pixel,  which  can  be 
analyzed  or  directly  displayed  in  tri-color  as  shown  in  Figure  17. 

By  employing  domain  decomposition  techniques,  the  coarse  resolution  subimages 
can  be  interpolated  and  summed  to  closely  approximate  a  fine  resolution  image  re¬ 
constructed  from  the  full  aperture,  fullband  SAR  signal.  Therefore,  phase  history 
decomposition  methods  are  computationally  efficient.  In  addition,  DAP  character¬ 
istics  of  bright  spots  or  peaks  can  be  automatically  extracted  from  the  subimages 
and  succinctly  presented  within  the  image  context  without  the  need  for  image  seg¬ 
mentation,  spectral  analysis,  or  most  importantly,  human  supervision.  Finally,  note 
that  polarimetric  analyses  commonly  performed  on  images  are  also  valid  for  use  with 
subimages.  Because  of  the  advantages  in  operational  and  computational  efficiency 
and  the  fact  that  canonical  scatterers  often  comprise  objects  of  interest,  this  disserta¬ 
tion  develops  a  new  SAR  imaging  and  scatterer  classification  theory  based  on  phase 
history  decomposition. 
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Figure  17.  SAR  image  highlighting  the  dispersive  characteristics  of  subimage  pixels. 
Subbands  centered  on  frequencies  8.8  GHz,  9.4  GHz,  and  10  GHz  are  coded  as  red, 
green,  and  blue  channels,  respectively  [46]. 


3.2  Proposed  Solution 

The  new  SAR  imaging  and  scatterer  classification  theory  centers  around  two  hy¬ 
potheses: 

•  it  is  possible  to  locate  and  classify  canonical  scatterers  by  observing  the  inten¬ 
sities  of  subimage  pixels,  and 

•  phase  history  decomposition  makes  this  approach  to  classification  highly  effi¬ 
cient. 

The  first  hypothesis  requires  development  of  a  new  model  to  predict  the  intensity  of 
subimage  peaks  due  to  canonical  scatterers.  Such  models  are  called  peak  models.  For 
SAR  signals,  spatial  resolution  is  inversely  proportional  to  spectral  bandwidth.  This 
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is  a  manifestation  of  the  Gabor  limit  [60].  Therefore,  phase  history  decomposition 
methods  essentially  operate  by  sacrificing  spatial  resolution  information  to  obtain 
spectral  bandwidth  information.  Typically,  the  spectral  diversity  of  the  subimages 
is  limited  in  order  to  maintain  reasonable  precision  for  scatterer  localization.  As  a 
result,  only  slowly  varying  pixel  intensities  can  be  accurately  measured  by  subimage 
analysis. 

All  canonical  scatterers  have  a  slowly  varying  amplitude  response  to  changes  in 
frequency,  in  accordance  with  Equation  (12).  However,  only  canonical  point  scat¬ 
terers,  defined  as  having  L  «  0  in  Equation  (13),  have  a  slowly  varying  ampli¬ 
tude  response  to  changes  in  azimuth.  Canonical  point  scatterers,  such  as  a  trihe¬ 
dral,  are  characterized  by  the  fact  that  their  peak  intensities  are  located  at  a  single 
location,  or  point,  in  the  image,  similar  to  the  peak  of  a  digitized  sine  function 
[42,  45,  57,  77,  122,  123,  131,  152,  153], 

In  contrast,  distributed  canonical  scatterers,  defined  as  having  L  >  0  in  Equation 
(13),  have  a  rapidly  varying,  sinc-like  amplitude  response  in  azimuth  [62,  122],  Dis¬ 
tributed  canonical  scatterers,  such  as  the  common  dihedral  with  a  fold-line  oriented 
parallel  to  the  imaging  plane,  are  characterized  by  the  fact  that  their  energy  usu¬ 
ally  spreads  across  multiple  pixels.  This  is  representative  of  a  digitized  rect  function 
with  a  region  of  support  over  multiple  samples  and  ripple  in  accordance  with  Gibbs 
phenomenon  [107].  Thus,  distributed  scatterers  often  appear  as  a  set  of  in-line  peaks 
of  approximately  the  same  amplitude.  For  example,  in  Figure  16,  Oval  C  shows  a 
likely  example  of  a  distributed  canonical  scatterer,  while  Ovals  A  and  B  show  likely 
examples  of  canonical  point  scatterers. 

Peak  models  have  already  been  developed  to  predict  the  intensity  of  subimage 
peaks  due  to  canonical  point  scatterers  [5,  57,  81,  117,  147,  150].  However,  no  model 
exists  to  predict  the  intensity  of  subimage  peaks  due  to  distributed  canonical  scatterers. 
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This  dissertation  develops  a  new  multi-peak  model  to  approximate  the  amplitudes  of 
localized  image  peaks  that  typically  appear  at  a  single  pixel  or  an  in-line  group  of  pix¬ 
els  in  a  SAR  image.  The  multi-peak  model  assumes  high  frequencies  and  uses  a  wide 
angle  approximation,  but  is  an  improvement  over  existing  peak  models  because  it 
explicitly  accounts  for  distributed  scatterers.  It  replaces  the  rapidly  varying  azimuth 
dependency  of  the  amplitude  function  in  the  spectral  domain  with  a  slowly  varying 
frequency  dependency.  In  this  way,  the  multi-peak  model  accounts  for  distributed 
scatterers  while  promoting  the  efficiencies  associated  with  phase  history  decomposi¬ 
tion  methods. 

The  multi-peak  model  for  a  canonical  scatterer  with  peak  amplitude  at  location 
(xq,  yq)  can  be  expressed  mathematically  as 

\Sqi.Xqi  Vq)  I  = 

where  wq  is  the  set  of  parameters  describing  the  amplitude  function  of  the  qth.  dis¬ 
tributed  canonical  scatterer,  while  w'q  is  a  reduced  set  of  parameters  describing  the 
amplitude  function  of  the  equivalent  canonical  point  scatterer.  The  equivalent  canon¬ 
ical  point  scatterer  has  an  azimuth-independent ,  scaled  amplitude  function  in  the 
spectral  domain  and  a  frequency  dependency  of  reduced  order.  The  approximation  is 
suitable  for  scatterers  of  sufficient  electrical  length  and  apertures  of  sufficient  width. 
These  are  the  high-frequency  and  wide-angle  assumptions  of  the  multi-peak  model. 
Chapter  IV  develops  the  multi-point  model  and  discusses  the  error  due  the  approxi¬ 
mation,  as  well  as  the  conditions  for  which  the  error  is  well-controlled. 

Not  all  image  peaks  are  due  to  canonical  scatterers.  Therefore,  there  must  be 
a  process  for  separating  image  peaks  due  to  canonical  scatterers  from  image  peaks 
due  to  non-canonical  scatterers.  Because  non-canonical  scatterers  have  amplitude 


B{HBHeSq(f,0;  w?)|  (xq,yq) 
B\HBHeSq(f;  w',)|  (x q,yq) 


(26) 
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responses  that  are  rapidly  varying  or  discontinuous  in  frequency,  their  image  peaks 
tend  not  to  persist  in  subimages  reconstructed  from  diverse  subbands.  Conversely, 
canonical  scatterers  have  a  slowly  varying  amplitude  response  in  frequency  and  tend 
to  persist.  Therefore,  the  SPLIT  algorithm  uses  a  persistence  criterion  to  reject  non- 
canonical  scatterers  from  the  scatterer  classification  process. 

For  each  persistent  peak,  SPLIT  compares  the  peak  intensities  from  all  subimages 
to  those  predicted  by  the  new  multi-peak  model.  For  example,  in  the  case  of  two 
subimages,  the  ratio  of  subimage  peak  intensities  is  related  to  the  subband  center 
frequencies  by  [57] 

\Mw*)\2  ^  (uy'+2  (27) 

\Mxq,vq)\2~\Uj  ’  1  J 

where  a'  is  the  parameter  describing  the  order  of  the  frequency  dependency.  Hence, 
the  value  of  ol  is  estimated  from  the  sampled  image  intensities  on  the  left-hand  side 
and  known  center  frequencies  on  the  right-hand  side.  Alternately,  when  more  than 
two  subimages  at  different  subband  center  frequencies  are  analyzed,  the  estimation  of 
a'  becomes  a  straightforward  curve  fitting  exercise  [57].  When  multiple  polarizations 
are  available,  the  SPLIT  algorithm  also  extracts  the  Krogager  parameters  from  the 
subimages.  A  least  squares  classifier  compares  the  extracted  parameters  to  the  ideal 
parameters  for  each  type  of  canonical  scatterer,  as  listed  in  Table  5.  The  SPLIT 
algorithm  and  the  error  due  to  the  narrow  band  approximation  in  Equation  (27)  are 
developed  in  detail  in  Chapter  V. 

The  multi-peak  model  and  SPLIT  algorithm  support  the  first  hypothesis,  while 
the  second  hypothesis  is  supported  by  the  integrated  algorithm.  The  integrated  al¬ 
gorithm  combines  SPLIT-based  classification  with  a  domain  decomposition  imaging 
algorithm.  The  result  is  a  computationally  efficient  scatterer  classification  algorithm 
that  autonomously  presents  classification  results  within  the  image  context.  The  com¬ 
putational  complexity  and  cost  of  the  integrated  algorithm  are  presented  in  detail  in 
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Table  5.  Extended  Frequency  Response  and  Polarimetric  Parameters  for  Ideal  Canon¬ 
ical  Scattering  Geometries. 


Scattering  geometry 

ex' 

K-q 

trihedral 

2 

1 

0 

dihedralgo 

2 

0 

l 

cylinder^ 

1 

l 

0 

top  hat 

1 

0 

1 

sphere,  plate 

0 

l 

0 

edge/wirego 

0 

0.5 

0.5 

dihedralo 

0 

0 

1 

cylinder0 

-1 

1 

0 

edge/wire0 

-2 

0.5 

0.5 

helical 

any 

0 

0 

Chapter  VI.  Furthermore,  all  SAR  imagery  for  the  experiments  in  Chapter  V  were 
produced  using  the  integrated  algorithm. 

In  summary,  the  multi-peak  model,  SPLIT  algorithm,  and  integrated  algorithm 
form  the  foundation  of  a  new  theory  for  efficiently  classifying  canonical  scatterers 
through  phase  history  decomposition.  The  operational  efficiency  is  the  key  moti¬ 
vation  where  DAP  characteristics  not  readily  available  in  the  SAR  image  are  made 
available  automatically  in  a  computationally  efficient  way.  These  characteristics  can 
be  succinctly  displayed  within  the  context  of  the  SAR  image  and  have  the  potential 
to  highlight  areas  of  interest  to  aid  in  the  management  of  high  precision  tools  and 
algorithms  for  SAR  image  analysis. 
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IV.  New  Multi-Peak  Model 


This  chapter  presents  the  new  multi-peak  model  to  approximate  the  amplitudes  of 
localized  image  peaks  that  typically  appear  at  a  single  pixel  location  or  as  an  in-line 
set  of  pixels  in  a  SAR  image.  The  multi-peak  model  is  derived  from  a  wide-angle 
approximation  of  the  well-known  parametric  model  for  canonical  scatterers  [62,  121], 
which  in  turn,  is  based  on  physical  models  of  electromagnetic  scattering.  In  this  way, 
the  multi-peak  model  is  an  improvement  over  existing  peak  models  which  poorly 
represent  distributed  canonical  scatterers,  such  as  the  common  dihedral  with  a  fold¬ 
line  oriented  parallel  to  the  imaging  plane. 

The  multi-peak  model  approximates  the  image  peak  amplitudes  due  to  distributed 
canonical  scatterers  as  if  they  are  due  to  an  equivalent  point  scatterer  with  an 
azimuth- independent,  dispersive  amplitude  function  in  the  spectral  domain.  The 
approximation  results  from  the  action  of  the  imaging  operator,  which  integrates  the 
sinc-like  reflectivity  pattern  in  azimuth  over  a  sufficient  aperture  width.  The  relative 
error  clue  to  the  approximation  is  shown  to  be  two  percent  or  less  when  a  tapered 
window  is  used  in  azimuth,  canonical  scatterers  are  ten  wavelengths  long  or  longer, 
and  aperture  widths  are  ten  degrees  wide  or  wider.  In  addition,  because  scatterers 
with  tilt  angles  near  0°  behave  as  distributed  scatterers  and  those  with  tilt  angles 
near  90°  behave  as  point  scatterers,  another  advantage  of  the  multi-peak  model  is 
that  scatterer  tilt  angles  near  0°  and  near  90°  can  be  discriminated  without  the  need 
for  fully-polarimetric  SAR  data. 

The  chapter  is  organized  as  follows.  Section  4.1  gives  the  mathematical  expres¬ 
sion  for  an  image  peak  due  to  a  canonical  scatterer.  It  discusses  how  and  under 
what  conditions  the  imaging  operator  integrates  the  sinc-like  amplitude  response  in 
azimuth  to  produce  a  point-like  amplitude  response  for  distributed  scatterers.  Sec¬ 
tion  4.2  presents  numerical  analysis  to  verify  that  error  clue  to  the  approximation 
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is  two  percent  or  less  when  canonical  scatterers  are  ten  wavelengths  long  or  longer 
and  aperture  widths  are  ten  degrees  wide  or  wider.  Section  4.3  presents  asymptotic 
analysis  which  further  reveals  that  the  model  error  is  best  controlled  with  the  use  of 
a  tapered  window  in  azimuth,  such  as  the  raised  cosine  windows  employed  in  SAR 
imaging  algorithms.  Section  4.4  presents  the  new  multi-peak  model  and  describes 
how  the  model  provides  a  reduced  feature  vector  for  scatterer  classification  by  phase 
history  decomposition.  Section  4.5  presents  some  additional  considerations  for  using 
the  multi-point  model  with  stripmap  mode  SAR  collection  geometries,  where  the  syn¬ 
thetic  aperture  is  usually  limited  to  a  few  degrees.  Finally,  Section  4.5.1  summarizes 
the  benefits  and  limitations  of  the  new  multi-peak  model. 

4.1  A  Peak  Model 

From  the  scatterer  models  in  Equations  (7)  of  Section  2.1.3,  (12)  of  Section  2.2.1, 
and  (13)  of  Section  2.2.2,  the  phase  history  due  to  a  canonical  scatterer  is  parame¬ 
terized  by 

w„)  =  S/(f;Aq,aq)  x  Se(f,0;L„0ai)e-™r'l  (28) 

where  Lq  >  0  for  a  distributed  canonical  scatterer  and  Lq  «  0  for  a  canonical  point 
scatterer,  setting  Sg  ~  1.  Thus,  from  Equations  (9)  and  (10)  of  Section  2.1.3,  the 
image  due  to  the  gth  canonical  scatterer  is 

•1  /»00  /»7T 

sq{x,y)  =  —  j  /  HB(f  -  fc)He(0  -  0c)Sf(f;Aq,aq) 

^  J  —  OO  J  —  7T  (29) 

xSg(f,e;Lq:e0q)e^ 
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The  peak  amplitude  of  Equation  (29)  resides  at  location  (xq,  yq)  in  the  image  where 
|r|  =  |r?|,V$.  Thus,  the  peak  amplitude  is 


Sq(Xq,  Vq)  ~  ^ 


HB(f-fc)Sf(f;Aq,aq)\f\ 

OO  _ 

x  f'  He(9  -  9-  L„90,)d9df. 


(30) 


Note  that  the  integral  in  azimuth  of  Equation  (29)  represents  an  inverse  Fourier 
transform  which  affects  the  image  of  point  scatterers  and  distributed  scatterers  dif¬ 
ferently.  A  point  scatterer  with  L  «  0  has  a  constant  or  slowly  varying  amplitude 
response  in  azimuth,  and  the  image  will  contain  a  single  peak  in  accordance  with 
Equation  (30).  However,  a  distributed  scatterer  with  L  0  has  a  sinc-like  amplitude 
function  in  azimuth,  and  given  a  sufficient  aperture  width,  the  image  will  contain  a 
rectangular  function  with  ripple.  The  ripple  is  a  manifestation  of  Gibbs  phenomenon 
[107]  and  results  in  one  or  more  image  peaks,  depending  upon  the  length  of  the  scat¬ 
terer  and  the  dimensions  of  the  image  pixels.  In  this  case,  Equation  (30)  predicts 
the  amplitude  of  the  rectangular  function,  including  any  peaks  caused  by  ripple.  The 
multi-peak  model  derives  its  namesake  from  this  effect,  where  distributed  canonical 
scatterers  typically  appear  as  an  in-line  group  of  multiple  peaks  in  the  image. 

Consider  a  distributed  canonical  scatterer  oriented  so  that  the  main  lobe  of  the 
sine  function  is  contained  within  the  azimuth  window.  In  this  case,  the  resultant  of  the 
inner  integral  of  Equation  (30)  is  dominated  by  the  area  under  the  main  lobe.  Indeed, 
the  sidelobes  are  small,  diminished  by  the  azimuth  window,  and  add  destructively,  so 
as  to  contribute  very  little  to  the  result.  Note  also  that  the  width  of  the  main  lobe 
is  inversely  proportional  to  Thus,  for  a  distributed  canonical  scatterer  of  fixed 
physical  length  and  oriented  so  that  the  main  lobe  of  the  sine  function  is  contained 
within  the  azimuth  window,  the  resultant  of  the  inner  integral  has  a  magnitude  that 
is  inversely  proportional  to  frequency. 
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This  idea  can  be  expressed  by  the  mathematical  approximation 


Hq(9  —  9c)Sg(f,  6]  L,  90)dO 


cA 

tv  >  yJmint 


(31) 


where  Omm  is  the  minimum  aperture  width  for  which  the  area  under  the  main  lobe  of 
the  sine  function  dominates.  Note  that  the  A  dependency  on  He  and  90  is  suppressed 
in  the  notation.  This  is  a  matter  of  convenience  due  to  the  assumption  that  azimuth 
window  and  orientation  angle  are  fixed  during  phase  history  decomposition. 

Note  that  the  approximation  in  Equation  (31)  depends  upon  the  electrical  length 
of  the  scatterer,  but  this  parameter  is  not  usually  controlled  by  the  radar  engineer. 
Thus,  an  appropriate  minimum  of  L  >  10A  is  chosen  in  accordance  with  the  high- 
frequency  approximations  of  the  underlying  GO/GTD  models  [83].  Note  also  that 
the  approximation  depends  upon  aperture  width  to  capture  the  mainlobe  of  the  sine 
function,  and  azimuth  window  to  suppress  the  sidelobes  of  the  sine  function.  Both 
of  these  are  controllable  by  the  radar  engineer  and  are  of  primary  importance.  The 
error  introduced  by  the  approximation  and  the  window  and  aperture  conditions  for 
which  the  error  is  well-controlled  are  best  understood  and  illustrated  by  numerical 
and  asymptotic  analysis. 


4.2  Numerical  Analysis 

Numerical  integration  of  Equation  (31)  for  a  range  of  frequencies  and  aperture 
widths  is  shown  in  Fig.  18,  for  the  case  of  L  —  1  meter  and  a  Hanning  window  for 
Hq.  Figures  18(a)  and  (b)  reveal  that  A  becomes  independent  of  frequency  when 
0  >  10°  and  /  >  3  GHz  for  this  case.  Figure  18(c)  also  shows  that  this  criteria 
is  met  for  varying  orientation  angles  6q.  Because  raised  cosine  windows,  other  than 
Hanning,  are  often  used  in  SAR  imaging  [76],  these  were  also  investigated  and  were 
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Figure  18.  Numerical  results  for  Eq.  (31)  where  H@  is  a  Hanning  function  with  region 
of  support  equal  to  0.  The  distributed  scatterer  is  of  length  L  =  1-meter. 


observed  to  produce  similar  results  to  that  shown  in  Fig.  18.  The  accuracy  of  the 
approximation  decreases  gradually  as  the  scatterer  becomes  electrically  shorter,  that 
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is  as  /  decreases  or  as  L  decreases  in  units  of  wavelength.  The  wavelength  is  defined 
as  A  —  c/  f. 

Because  the  model  in  Equation  (29)  is  based  on  GO/GTD,  it  is  most  accurate 
for  electrically  large  scatterers,  where  L  >  10 A.  This  is  sometimes  referred  to  as  the 
physical  optics  region  or  high  frequency  approximation  for  electromagnetic  scattering. 
Accordingly,  the  numerical  integration  results  in  Figure  18  are  truncated  at  fmin  =  3 
GHz,  where  the  L  —  1  meter  scatterer  is  ten  wavelengths  long.  Therefore,  (~)min  ~  10° 
is  sufficient  to  support  a  high  frequency  approximation  for  electromagnetic  scattering, 
and  for  the  purposes  of  this  paper,  apertures  greater  than  ten  degrees  are  considered 
to  constitute  wide-angle  SAR.  In  other  contexts,  the  definition  for  wide-angle  SAR 
is  constrained  to  the  condition  where  bandwidth  in  frequency  is  much  more  limited 
than  bandwidth  in  azimuth  [96,  110,  146],  but  Equation  (31)  is  not  limited  to  this 
condition. 


4.3  Asymptotic  Analysis 


Setting  the  limits  of  integration  in  Equation  (31)  to  the  region  of  support  for  Hq, 
9  e  [9C  —  ® ,  9C  +  ® ],  with  the  constraint  ®  <  |  —  |0O|  and  setting  9C  =  0  without  loss 
of  generality  gives 

He(9)  sine  [|  fL  8111(6*  —  d0)]  d,9.  (32) 

The  change  of  variables  a  =  and  x  =  sin(d  —  9 o)  further  simplifies  the  integral  to 


A 


i-i 


Hq(90  +  sin  1  x) 


\/l  — 


x- 


a  sine  (ax)dx, 


(33) 


where  the  earlier  constraint  causes  the  limits  of  integration  to  be  constrained  to  the 
interval  [sin  (—  ®  —  90)  ,  sin  (®  —  6*0)] ,  or  more  generally  x  e  [—1, 1],  As  a  approaches 
infinity,  the  second  factor  in  Equation  (33)  is  a  form  of  the  Dirac  delta  function  [95]. 


56 


Under  this  limit,  the  sine  function  becomes  a  sampling  function  so  that  the  first  factor 
is  sampled  at  x  —  0,  and 


lim  A 

a— kx) 


Hq(90  +  sin  1  x) 
a/1  —  x 2 


a  sine (ax)dx 


=  He(90), 


(34) 


where  Hq(90)  is  assumed  positive  for  |0O|  <  §  and  zero  otherwise.  Under  the  limit, 
A  equals  a  constant;  therefore,  the  error  introduced  by  the  wide-angle  approximation 
in  Equation  (31)  is  bounded.  That  is,  for  any  typical  azimuth  window,  the  error 
bound  diminishes  with  increasing  frequency,  or  equivalently  as  the  distributed  scatter 
increases  in  electrical  length. 

Bounded  error  alone  is  not  sufficient  to  make  the  wide-angle  approximation  gen¬ 
erally  useful  for  time-frequency  analysis  of  SAR  imagery.  It  is  also  important  that  A 
be  insensitive  to  changes  in  frequency,  or  equivalently,  changes  in  a.  This  sensitivity 
is  revealed  by  taking  the  partial  derivative  with  respect  to  a  of  the  right  hand  side  of 
Equation  (33),  where  a  partial  derivative  equal  to  zero  reveals  that  A  is  independent 
of  a.  Thus,  noting  that  9('abnA"J  11  —  cos(7raa;),  the  wide-angle  approximation  relative 
error  is  defined  as 


e(He,a,90 ) 


Hstfo+sin  1  x) 
\/l—x2 


cos(nax)dx 


He(90) 


(35) 


where  e  =  0  indicates  that  A  is  independent  of  a.  Plots  of  the  error  are  shown  in 
Fig.  19  for  the  case  of  90  =  0  and  for  cases  a  =  20  and  a  =  40.  The  plots  compare 
the  errors  of  a  rectangular  azimuth  window  and  three  other  windows  commonly  used 
in  SAR  image  processing  [76].  All  windows  are  in  accordance  to  Matlab®  default 
definitions.  The  Hanning  window  appears  to  perform  best  overall,  while  the  Taylor 
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0  (degrees)  0  (degrees) 

Figure  19.  Wide-Angle  approximation  relative  error  for  a  rectangular  window  in  az¬ 
imuth  and  three  other  windows  commonly  used  in  SAR  image  processing. 


window  performance  is  marginal,  and  the  rectangular  window  performance  is  notably 
poor.  Similar  results  are  obtained  for  larger  values  of  a,  but  with  lower  overall  error, 
as  expected.  Also  as  expected,  the  error  increases  rapidly  for  an  aperture  width  less 
than  ten  degrees,  regardless  of  window  type.  In  addition,  the  results  are  similar  for 
varying  values  of  6q.  Because  a  rectangular  window  in  azimuth  causes  much  larger 
error  than  that  caused  by  a  tapered  window,  the  multi-peak  model  is  recommended 
for  use  with  tapered  windows,  such  as  the  family  of  raised  cosine  windows  typically 
employed  for  SAR  imagery. 


4.4  The  Multi-Peak  Model 

Based  on  the  azimuth  limitations  of  the  model,  the  peak  amplitude  of  the  qth 
scattering  center  in  Equation  (30)  is  approximated  as 


sq{xq,yq)  «  —  /  HB(f  ~  fc)Aq(jf)a<l/2\f\ 


cAq 
2  fLq 


df 


(36) 


58 


for  scatterers  of  sufficient  length  (L  >  10A)  and  an  aperture  of  sufficient  width  (0  > 
10°).  Here,  the  phase  history  due  to  the  qth  scatterer  is  independent  of  6.  This  result 
is  derived  from  Equation  (30)  when  expressed  in  polar  coordinates,  but  a  similar 
result  can  be  obtained  for  rectangular  coordinates  using  the  far  field  approximation, 
R(x,y )  ~  xcosd  +  ysmO,  and  the  small  angle  approximation,  /  =  >J +  fy  ~ 
fy.  Therefore,  for  rectangularly  formatted  phase  histories,  an  additional  aperture 
limitation  of  0  <  20°  is  required  in  order  to  support  a  small  angle  approximation. 

Because  L  is  usually  not  known  a  priori ,  it  is  desirable  to  modify  Equation  (36) 
to  account  for  both  distributed  and  point  scatterers  with  a  single  parameter,  a'. 
Adopting  the  form  of  Equation  (12),  a  convenient  form  for  Equation  (36)  is 


Sq{Xq>  Dq)  ~  ^ 


HbU  -  fc)He(0  -  0c)  kuf)“i/2l  I/I m,  (37) 


—  OO  J  —  7T 


where 

Sq(f-,W,)  =  A'q(jf)</2e-*k^  (38) 

is  the  phase  history  for  an  equivalent  point  scatterer  to  Equation  (28)  of  Section  4.1. 

For  canonical  point  scatterers,  L  zn  0,  which  in  turn  causes  a'  —  a  and  A'  =  A.  For 
distributed  canonical  scatters,  L  >  10A  and  0  >  10°,  which  in  turn  causes  a'  =  a  —  2 
and  A'  =  Ho(0-0  )de-  This  supports  the  hypothesis  expressed  in  Equation  (26) 

of  Section  3.2  where  w  =  [x,  y,  A,  a,  L,  0o\  from  Equation  (28)  and  w;  =  [x,y,A!,a!] 
from  Equation  (38).  Therefore,  the  SAR  image  peak  amplitude  clue  to  a  distributed 
canonical  scatterer  can  be  modeled  as  clue  to  an  equivalent  canonical  point  scatterer 
having  an  azimuth- independent,  scaled  amplitude  function  in  the  spectral  domain 
and  a  frequency  dependency  of  reduced  order. 

The  fact  that  image  peaks  due  to  any  canonical  scatterer  can  now  be  approximated 
as  clue  to  a  canonical  point  scatterer  is  central  to  the  high-frequency  multi-peak  model 
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Table  6.  Extended  Frequency  Response  Parameter  for  Ideal  Canonical  Scattering  Ge¬ 
ometries. 


Scattering  geometry 

ex' 

trihedral,  dihedral90 

2 

cylinder 90,  top  hat 

1 

sphere,  plate,  edge/wire90,  dihedral0 

0 

cylinder0 

-1 

edge/wireo 

-2 

for  wide-angle  SAR.  This  new  model  provides  an  extension  of  the  traditional  point 
model  in  Equation  (12)  of  Section  2.2.1,  and  as  a  result,  Table  2  of  Section  2.2.1 
expands  to  become  Table  6.  Aside  from  the  plate,  which  is  symmetric  under  the 
model  in  two  dimensions,  the  distributed  scatterers  are  discriminated  from  their  point 
scatterer  variants  by  tilt  angles  of  0°  and  90°,  as  indicated  by  a  subscript.  A  benefit 
of  the  multi-peak  model  is  that  scatterer  tilt  angles  of  r  «  0°  and  r  ~  90°  can  be 
discriminated  without  the  need  for  fully-polarimctric  SAR  data. 

If  fully  polarimetric  data  is  available,  ambiguity  between  most  scatters  can  be 
further  resolved  using  a  Krogager  decomposition  to  obtain  proportions  of  the  odd- 
bounce,  even-bounce,  and  helical  scattering  [86].  Table  5  lists  the  new  frequency 
parameter  a'  with  the  odd  and  even  bounce  parameters  kq  and  ne,  respectively.  Us¬ 
ing  this  basis  for  classification,  only  the  sphere  and  plate  are  ambiguous.  When  the 
phase  history  is  imaged  over  multiple,  wide-angle  subapertures,  then  for  reasonable 
elevation  angles  (e.g.  less  than  60°),  point  scatterers  will  tend  to  persist  across  mul¬ 
tiple  subaperture  images,  whereas,  distributed  scatterers  do  not  persist.  In  this  case, 
a  persistence  criteria  can  be  used  to  further  resolve  ambiguity  between  distributed 
and  point  scatterers,  such  as  the  sphere  and  plate.  The  papers  [26,  96,  110,  137]  and 
their  references  contain  information  on  scatterer  persistence  which,  for  brevity,  are 
not  addressed  in  this  dissertation. 
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4.5  Stripmap  Mode  SAR  Considerations 


The  multi-peak  model  is  accurate  for  the  case  of  L  >  10A  and  0  >  10°.  However, 
stripmap  mode  SAR  systems  often  have  aperture  widths  less  than  10°.  In  addi¬ 
tion,  the  multi-peak  model  also  assumes  an  imaging  operator  that  is  a  2D  transform 
between  the  spectral  and  spatial  domains.  However,  some  stripmap  mode  imaging  al¬ 
gorithms  do  not  operate  in  the  2D  spectral  domain.  In  response  to  these  restrictions, 
this  section  examines  the  special  considerations  that  govern  the  use  of  the  multi-peak 
model  with  stripmap  mode  SAR  systems. 

The  angular  diversity  in  stripmap  SAR  signals  is  ultimately  restricted  to  the  width 
of  the  azimuth  beamwidth.  Furthermore,  a  narrow  azimuth  beamwidth  is  usually 
desired  for  two  primary  reasons.  The  first  is  a  matter  of  antenna  gain.  A  narrow 
beamwidth  creates  a  higher  antenna  gain,  improving  the  overall  signal  to  noise  ratio 
of  the  SAR  system  [36].  The  second  is  a  matter  of  pulse  repetition  frequency  (PRF). 
The  PRF  determines  the  Doppler  sampling  rate,  and  so  a  higher  PRF  provides  for 
a  larger  angular  diversity  without  aliasing  [36].  However,  the  PRF  should  also  be 
low  enough  that  the  signal  corresponding  to  the  entire  swath  width  can  be  received 
without  range  ambiguity  [36].  These  two  requirements  are  in  competition,  and  the 
PRF  is  usually  chosen  to  limit  angular  diversity  so  as  to  obtain  a  wider  swath  width. 
Thus,  as  angular  diversity  is  limited  by  PRF,  very  large  azimuth  beamwidths  provide 
no  additional  advantage,  in  most  cases. 

The  half-power  beamwidth  of  a  radar  system  can  be  approximated  by  [36] 

0  ~  0.88-^j,  (39) 

JcL 

where  0  is  in  units  of  radians,  c  is  the  speed  of  light,  fc  is  the  center  frequency  of 
the  radar  signal,  and  l  is  the  physical  length  of  the  antenna.  For  convenience,  typical 
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Table  7.  Azimuth  Beamwidths  of  Various  SAR  System  Types  [52,  36,  37,  148,  129,  98, 
128,  132,  97,  114,  69,  143]. 


© 

System  Type 

Band 

l(  m) 

<  1° 

spaceborne 

X 

5-15 

<  1° 

spaceborne 

C 

8-15 

1° 

spaceborne 

L 

8-15 

1°  -  2° 

airborne 

X 

0.6-1. 1 

2°  -  4° 

airborne 

C 

0. 7-1.4 

7°  -  15° 

airborne 

L 

0.7-1. 6 

20°  -  115° 

airborne 

P 

0.7-1. 8 

azimuth  beamwidths  for  a  variety  of  SAR  systems  are  listed  in  Table  7.  Here, 
typical  wavelengths  for  X-,  C-,  L-band  radars  were  approximately  9.6,  5.3,  1.3  GHz, 
respectively,  while  P-band  frequencies  varied. 

For  stripmap  mode  SAR,  only  airborne  systems  in  L-band  and  P-band  are  ex¬ 
pected  to  provide  a  large  angular  diversity,  while  other  systems  typically  have  an 
angular  diversity  of  only  a  few  degrees  or  less.  For  convenience,  the  remaining  dis¬ 
cussion  equates  the  angular  diversity  or  aperture  of  the  SAR  system  to  its  half-power 
beamwidth  in  azimuth. 

One  significant  limitation  of  the  multi-point  model  is  that  it  assumes  azimuth 
integration  is  performed  in  the  two-dimensional  (2D)  spectral  domain.  While  in¬ 
tegration  of  2D  spectral  domain  data  is  standard  for  most  spotlight  mode  imaging 
algorithms,  it  is  not  standard  for  stripmap  mode  imaging  algorithms.  For  instance, 
the  range-Doppler  algorithm  (RDA)  and  the  chirp  scaling  algorithm  (CSA)  perform 
azimuth  integration  in  the  range-Doppler  domain,  not  the  2D  spectral  domain  [36]. 
The  range-Doppler  domain  is  obtained  by  performing  range-frequency  integration 
before  azimuth  integration.  As  a  result,  the  azimuth  response  of  canonical  scatter¬ 
ed  is  not  integrated  out  at  an  appropriate  level  to  allow  the  amplitude  response  of 
distributed  canonical  scatterers  to  become  azimuth-independent.  Therefore,  a  key 
assumption  of  the  new  multi-point  model  is  violated  for  these  imaging  algorithms. 
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Figure  20.  Relative  wide-angle  approximation  errors  for  distributed  canonical  scatter- 
ers  of  varying  lengths.  Hq  is  a  Kaiser  window  with  a  coefficient  equal  to  3.6  and  a 
region  of  support  equal  to  0. 


Fortunately,  the  omega-K  algorithm  (OKA)  is  a  stripmap  mode  imaging  algorithm 
that  has  the  advantage  of  processing  data  in  the  2D  spectral  domain.  In  addition, 
the  OKA  is  considered  the  most  accurate  of  the  stripmap  mode  imaging  algorithms 
for  SAR  data  collected  over  a  wide-aperture  [36].  For  these  reasons,  the  OKA  is 
recommended  for  use  with  the  new  scatterer  classification  method  for  wide-angle 
SAR,  whereas  the  RDA  and  OS  A  are  not. 

The  accuracy  of  the  multi-point  model  degrades  with  a  narrowing  aperture  or 
a  reduction  in  the  electrical  length  of  the  scatterer.  Specifically,  as  the  scatterer 
becomes  electrically  longer,  the  model  tolerates  a  narrower  aperture.  For  example, 
plots  of  the  relative  wide-angle  approximation  error  from  Equation  (35)  for  a  Kaiser 
window  with  a  coefficient  equal  to  3.6,  an  orientation  angle  of  zero  degrees,  and  a 
range  of  scatterer  lengths  are  provided  in  Figure  20.  (The  Kaiser  window  simulates 
the  azimuth  beam  pattern  [36].)  It  can  be  seen  that  for  a  given  subaperture  width, 
the  minimum  length  for  which  the  multi-point  model  error  is  kept  at  near  two-percent 
or  less  is  Lmin  &  8OA/0,  where  0  is  expressed  in  degrees.  Using  this  error  threshold 
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as  an  example  and  combining  these  results  with  Eq.  39,  produces  the  relationship 


Lmin  «  1.6 1,  (40) 

which  can  be  used  as  a  rule  of  thumb  for  estimating  the  minimum  effective  length 
of  distributed  scatterers  which  will  be  well  modeled  for  a  given  SAR  system.  Addi¬ 
tional  analysis  reveals  that  the  factor  1.6  in  Equation  (40)  is  directly  related  to  the 
broadening  factor  of  the  azimuth  window.  It  follows  that  a  window  with  a  smaller 
broadening  factor  will  produce  a  first  null  in  the  error  function  at  a  smaller  aperture 
width,  but  at  the  expense  of  higher  sidelobes.  Sidclobes  in  the  error  function  are 
controlled  by  choice  of  azimuth  window.  Additional  analysis  revealed  that  the  raised 
cosine  windows  typically  employed  in  radar  imaging,  such  as  Taylor,  Hanning  and 
Hamming,  are  sufficient  to  control  modeling  error.  Such  windows  all  have  similar 
broadening  factors  which  are  not  expected  to  have  a  great  effect  on  the  relationship 
given  in  Equation  (40). 

Using  Equation  (40)  to  evaluate  the  SAR  systems  listed  in  Table  7  reveals  that 
the  spaceborne  radars  possess  antenna  lengths  which  require  distributed  canonical 
scatterers  with  effective  lengths  of  10  to  20  meters  in  order  to  ensure  accuracy  of  the 
multi-point  model.  In  contrast,  airborne  radars  are  an  order  of  magnitude  shorter 
and  effective  lengths  of  only  a  few  meters  are  sufficient.  Last,  note  that  the  above 
analysis  focused  on  distributed  canonical  scatterers,  with  an  effective  length  L  >  0. 
Recall  that,  canonical  point  scatterers,  such  as  spheres  and  trihedrals,  have  a  constant 
or  slowly  varying  angular  response  in  azimuth.  Therefore,  canonical  point  scatterers 
have  an  effective  length  of  L  «  0  and  are  modeled  accurately,  regardless  of  aperture 
width  or  antenna  size. 
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4.5.1  Summary  of  Benefits  and  Limitations. 


This  chapter  presented  a  new  multi-peak  model,  as  expressed  in  Equations  (37) 
and  (38)  of  Section  4.4,  to  predict  the  dispersive  behavior  of  subimage  peaks  due 
to  canonical  scatterers.  The  model  was  shown  to  be  an  improvement  over  existing 
peak  models  because  it  accounts  for  common  distributed  canonical  scatterers,  includ¬ 
ing  plates  at  broadside  aspect,  as  well  as  dihedrals,  cylinders,  and  edges/wires  lying 
parallel  to  the  imaging  plane.  For  distributed  canonical  scatterers,  the  model  ap¬ 
proximates  the  amplitude  response  with  a  single  factor  having  an  inverse  frequency 
dependency.  This  factor  replaces  the  integration  of  the  mainlobe  of  the  sinc-like  re¬ 
flectivity  pattern  in  azimuth  as  part  of  the  imaging  process.  The  approximation  was 
shown  to  be  valid  for  wide-angle  SAR  when  the  synthetic  aperture  is  greater  than  ten 
degrees  in  azimuth.  Furthermore,  the  model  accuracy  was  shown  to  degrade  gradually 
with  a  decrease  in  aperture  width  below  ten  degrees.  The  wide-angle  approximation 
error  was  shown  to  be  controllable,  but  is  quite  high  for  a  rectangular  window  in 
azimuth.  Therefore,  the  new  multi-peak  model  is  recommended  for  use  with  tapered 
windows. 

The  new  multi-peak  model  is  useful  for  scatterer  classification  by  phase  history 
domain  decomposition  methods.  These  methods  sacrifice  some  precision  in  scatterer 
localization  in  order  to  gain  efficiency  by  removing  the  need  for  human  supervision 
typically  required  for  image  segmentation  methods.  Future  follow-on  research  may 
reveal  even  more  benefits  of  the  new  multi-peak  model.  However,  notable  limitations 
of  the  model  include  the  following.  Accuracy  of  the  model  is  limited  to  the  restricted 
apertures  0  >  10°  for  polar  formatted  data  and  10°  <  0  <  20°  for  rectangular 
formatted  data;  scatterers  must  be  of  sufficient  length  ( L  >  10A)  to  support  a  high- 
frequency  approximation;  and  a  tapered  window  must  be  employed  in  azimuth  to 
suppress  the  sidelobes  of  the  sine  function.  In  addition,  the  multi-peak  model  is  de- 
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veloped  from  GO/GTD  scattering  models  which  assume  a  perfect  electrical  conductor 
in  the  far  field.  Furthermore,  the  model  assumes  a  normalized  phase  history  account¬ 
ing  for  antenna  gain  pattern,  spherical  wave  propagation,  and  other  losses.  Lastly, 
for  stripmap  mode  SAR,  the  accuracy  of  the  model  is  limited  to  use  of  the  Ornega- 
K  algorithm  and  distributed  scatterers  having  an  effective  length  approximately  1.6 
times  the  physical  length  of  the  antenna. 
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V.  SPLIT  Algorithm 


This  chapter  develops  the  SPLIT  algorithm  for  classifying  canonical  scatterers. 
Section  5.1  explains  how  the  SPLIT  algorithm  analyzes  the  pixels  of  all  subimages 
within  a  given  subaperture  and  extracts  a  feature  vector  based  on  the  multi-peak 
model.  Section  5.2  describes  how  these  feature  vectors  are  summed  across  multi¬ 
ple  subapertures.  Section  5.3  develops  a  least  squares  classification  algorithm  using 
the  Euclidean  norm  for  a  set  of  ideal  feature  vectors,  each  representing  a  class  of 
canonical  scatterer.  Section  5.4  presents  SPLIT  scatterer  classification  results  for 
simulated  and  measured  data.  Finally,  Section  5.5  examines  the  sensitivity  of  SPLIT 
classification  accuracy  to  signal  parameters,  such  bandwidth  and  interference  due  to 
neighboring  canonical  scatterers,  clutter,  and  noise.  It  presents  a  statistical  signal 
model  with  associated  Monte  Carlo  simulations  to  illustrate  the  impact  of  bandwidth 
and  interference  on  classification  accuracy. 

5.1  Subaperture  Feature  Extraction 

The  SPLIT  algorithm  analyzes  the  pixels  of  all  subimages  within  a  given  sub- 
aperture  and  extracts  a  feature  vector  based  on  the  multi-peak  model.  Given  an 
M  x  N  subimage,  the  feature  vector  for  a  peak  at  the  subimage  pixel  (m,  n)  in  the 
j th  subaperture  is 

Vj(m,  n)  =  [a',  k0,  ne]j,  (41) 

where  a'  is  the  frequency  parameter  of  the  multi-peak  model  and  k0  and  ne  are  the 
proportion  of  odd-bounce  and  even-bounce  scattering  energy,  respectively. 

The  following  subsections  describe  the  steps  of  the  SPLIT  algorithm.  First,  the 
scattering  centers  are  localized  to  subimage  pixels  ( mq,nq )  using  a  peak  detection 
algorithm.  The  peak  detection  algorithm  identifies  any  peaks  appearing  at  the  same 
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Figure  21.  Scattering  centers  are  stationary  in  range  with  changes  in  frequency  [28]. 


location  in  all  subimages.  These  stable  peaks  indicate  a  stationary  and  slowly  varying 
amplitude  response  to  changes  in  frequency,  and  thus,  are  likely  due  to  canonical 
point  scatterers.  Next,  SPLIT  extracts  the  frequency  parameter  from  the  available 
co-polarization  channels  using  a  least  squares  estimator.  Finally,  SPLIT  extracts 
the  polarimetric  parameters  from  each  subimage  using  Krogager  decomposition,  as 
applicable. 

5.1.1  Scattering  Center  Localization. 

Scattering  centers  are  stationary  in  range  with  changes  in  frequency,  as  illustrated 
by  the  simulated  ID  range  profiles  of  Figure  21.  In  addition,  the  multi-peak  model 
predicts  that  image  peaks  due  to  canonical  point  scatterers  have  an  amplitude  that 
varies  slowly  with  changes  in  frequency.  Therefore,  the  subimage  peaks  due  to  canon¬ 
ical  point  scatterers  are  expected  to  be  stable  and  occur  at  the  same  pixel  coordinates 
in  each  subimage.  In  contrast,  non-canonical  scatterers,  such  as  those  due  to  reso¬ 
nance,  material  dispersion,  structural  dispersion,  or  scatter  motion,  are  not  expected 
to  produce  stable  subimage  peaks.  In  this  way,  the  SPLIT  algorithm  determines  the 
likely  locations  of  canonical  scatters  using  a  stable  subimage  peaks  detector. 

The  stable  subimage  peaks  detector  works  as  follows.  For  subimages  of  the  jth 
subaperture,  the  SPLIT  algorithm  assigns  a  ‘1’  to  each  pixel  that  is  a  local  peak 
as  compared  to  its  neighboring  pixels.  It  assigns  all  other  pixels  a  ‘0.’  The  SPLIT 


algorithm  arranges  the  pixel  values  in  an  M  x  N  array,  one  for  each  subimage.  An 
illustrative  example  featuring  three  subimages  is  given  by 
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where  the  output  of  the  peak  detector,  Z>{-},  for  each  subimage  is 


®{l9.j|2} 
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0  0  10 

10  0  0 

0  0  0  1 
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0  0  0  0 
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The  peak  detector  can  be  readily  implemented  using  the  imregionalmax  function 
in  Matlab®.  The  stable  peaks  are  found  by  a  simple  element-wise  multiplication, 
sometimes  referred  to  as  an  array  multiply.  The  result  of  the  stable  peak  detector  is 


0  0  0  0 


i 

0®{fei2} 


0  0  10 
10  0  0 


0  0  0  0 


(42) 


where  (£)  represents  array  multiplication  and  I  =  3  is  the  number  of  subimages.  For 
the  illustrative  example  in  Equation  (42),  there  are  only  two  stable  subimage  peaks, 
and  these  pixel  coordinates  reveal  the  likely  locations  of  canonical  scatterers. 
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5.1.2  Frequency  Parameter  Estimation. 


The  SPLIT  algorithm  obtains  a  frequency  parameter  estimate  by  examining  the 
intensities  of  the  stable  subimage  pixels  and  comparing  these  intensities  to  those  pre¬ 
dicted  by  the  multi-peak  model.  This  section  presents  a  mathematical  development 
for  a  simple  exponential  relationship  between  the  measured  subimage  pixel  intensities 
and  the  known  subband  center  frequencies.  The  exponent  is  directly  related  to  the 
frequency  parameter,  a'  of  the  multi-peak  model.  The  relationship  is  based  on  a  first 
order  approximation  using  a  narrow-band  assumption.  The  following  subsections  ex¬ 
amine  the  error  due  to  the  narrow-band  approximation  and  describe  a  curve  fitting 
procedure  for  estimating  the  frequency  parameter  from  the  peak  intensity  measure¬ 
ments. 


5. 1.2.1  First  Order  Approximation. 

Starting  with  the  multi-peak  model  for  the  qth  scattering  center  from  Equation 
(37)  of  Section  4.4  and  signifying  the  wide-angle  approximation  with  the  ^  symbol 
gives  the  equality 


^q{,XqiVq)  ^ 


He(9  -  9c)d9 


900  A/ 

HbW*  fc)-(jf)1+a',2df.  (43) 

-00  7 


Setting  C  =  f  He(9  —  9c)d9 ,  assuming  only  positive  frequencies,  and  applying  a 
Taylor  series  expansion  of  (j f)a'^2  in  the  neighborhood  of  fc  gives 


sq(xq,yq)=C  HR(f  -  fc) 


m™*  +  E  sr rUf)M'2 


n= 1 


dfr 


(/  -  fc y 


f=fc 


df. 

(44) 
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Defining  hB{t)  as  the  inverse  Fourier  transform  (IFT)  of  HB(f)  gives 


sq(xq,yq)  =  C(jfc)1+a'/2hB(  0)  + 

n— 1 


1  dn(jf)1+a '/2 
n\  dfn 


(45) 


where  the  second  term  contains  higher-order  functions  resulting  from  the  Taylor  series 
expansion  and  frequency  domain  differentiation  property  of  Fourier  transforms  [107]. 

Note  that  a  symmetric  azimuth  window  will  cause  all  odd-ordered  derivatives  of 
h^\ 0)  to  become  zero.  Furthermore,  for  a'  G  {—2,  0}  the  derivative  with  respect  to  / 
becomes  zero  for  all  n.  These  facts  combined  with  a  narrow-band  assumption,  causes 
the  summation  term  in  Equation  (45)  to  be  negligibly  small  in  practice.  For  example, 
in  the  case  of  4.0  GHz  bandwidth  at  X-band,  the  relative  error  due  to  the  approxi¬ 
mation  is  less  than  one-percent.  An  analysis  of  the  narrow-band  approximation  error 
is  provided  in  the  following  subsection. 

Consider  the  stable  subimage  pixel  ( mq,  nq )  dominated  by  the  response  due  to  a 
canonical  scatterer  at  coordinates  ( xq,yq ).  In  this  case,  a  first-order  approximation  of 
the  ratio  of  the  pixel  intensities  between  two  subimages,  having  subwindows  centered 
at  two  different  center  frequencies,  fc\  and  fc 2  is 


\gij(mq,nq)\2  \C(jfcl)1+a'/2hB{0)\2 
\g2j(rnq,nq)\2  ~  \C{j  fc2)l+a' !2hB{Q)\2' 


Assuming  a  uniformity  in  the  sampling  and  weighting  of  each  subdomain  (i.e.  sub¬ 
windows  in  azimuth  are  identical  and  each  subband  has  an  equal  number  of  uniform 
samples)  allows  the  constants  in  the  numerator  and  denominator  to  cancel.  This 
results  in  a  simple  relationship  between  the  measured  subimage  pixel  intensities  and 
the  subband  center  frequencies  given  by 


1 9ij(mq,nq)\2  _  / fci\a'+2 
\g2j(mq,nq)\2  \fc2J 


(47) 
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where  the  accuracy  of  the  approximation  is  subject  to  three  conditions: 

1.  the  energy  in  the  pixel  is  dominated  by  the  backscattering  from  a  canonical 
scatterer; 

2.  the  wide-angle  approximation  error  of  the  multi-peak  model  is  well-controlled; 
and 

3.  the  narrow-band  approximation  error  resulting  from  ignoring  the  higher-order 
terms  in  Equation  (45)  is  negligible. 

The  first  assumption  is  addressed  later  in  Section  5.5,  the  wide-angle  approxima¬ 
tion  error  was  previously  addressed  in  Sections  4.2  and  4.3,  and  the  narrow-band 
approximation  error  is  discussed  in  the  following  subsection. 


5. 1.2. 2  Narrow-band  Approximation  Error. 

This  section  examines  the  narrow-band  approximation  error,  which  results  from 
ignoring  the  higher-order  terms  in  Equation  (45).  The  following  derivation  assumes 
that  the  subwindow  in  frequency  is  symmetric,  with  frequencies  restricted  to  positive 
values,  so  that  fc>  0  and  B  <  2 fc.  Under  these  conditions,  the  band  limited  integral 
of  Equation  (43)  becomes 

rfc+B/ 2 

*,(*„  V,)=C  HB(f  -  fc)(jf)1+°'f2df.  (48) 

Jfc-B/2 


Next,  note  that 


(j/)I+“72  =  (v7)°+2/1+“72 


(49) 
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In  order  to  analyze  the  error  in  the  neighborhood  of  /c,  let 


(jf)1+a'/2  =  (jfc)1+a’/2  +  (jf)1+a'/2  ~  (jfc)1+a'/2 


(±0-  +3) 
\  V2 


(±(1  +  3) 

V  V2 


a'+2 

j'l+a' / 2  j\+a'/2  _  yl+a'/2 


a'+ 2 


jl+a' /2 


f 

1+1  J- 

Jc 


l+a'/2 


-  1 


(50) 


Choosing  to  restrict  the  function  (jf)1+a^2  to  positive  roots,  recalling  that  the  window 
is  centered  on  frequency  fc,  setting  D  =  C  fc+°  ^ and  substituting  into 

Equation  (48)  gives 


^  rfc+B/2 

sq(xq,yq)  =  D  HB(f  —  fc) 

Jfc-B/2 


f 

1+1  j- 

Jc 


l+a'/2 


-  1 


df 


(51) 


=  D  M  0)+  /  HB(f-fc ) 


'  fc—B/2 


l+a'/2 


df}, 


where  the  hrst  term  represents  the  Erst  order  approximation  of  Equation  (45)  used  in 
Equation  (46)  and  the  integral  represents  the  error  terms.  The  narrow-band  relative 
error  is  defined  as  the  magnitude  of  the  error  terms  divided  by  the  first  order  ap¬ 
proximation  term.  With  the  change  of  variables  u  =  f  —  fc,  the  narrow-band  relative 
error  is 


£(Hb,I3) 


1+0:72 


1  du 


hB(  0) 


(52) 


where  the  error  is  dependent  upon  the  window  function  and  fractional  bandwidth, 

P  =  B/U 
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A  bound  on  the  error  can  be  derived  by  examining  the  integral  in  the  numerator  of 
Equation  (52)  partitioned  about  zero.  Assuming  a  positive- valued  window,  HB(u )  > 
0,  then  for  the  case  of  a'  >  —2: 


Hb(u) 


'-B/2 


.  l+a'/2 

i  +  £l  -i 

Jc 


fB/ 2 

du  <  0  <  /  Hb(u ) 

Jo 


.  l+a'/2 

i  +  £l  -i 

Jc 


du 


(53) 

Alternately,  for  the  case  of  a'  <  —2,  the  inequalities  in  Equation  (53)  are  reversed. 
Note  also  that 


Hb(u) 


'-B/2 


,  U 

+  7c 


l+a'/2 


-  1 


du 


> 


cB/2 


Hb(u ) 


,  u 

+  Jc 


l+a'/2 


-  1 


du  . 
(54) 


Thus,  from  Equations  (51)  through  (54)  the  error  is  bounded  by 


e(HB,P)< 


-  sgn(a/  +  2) 

hB(  0) 


HB(u) 


'  —  B/2 


.  l+a'/2 

i  +  fi  -i 

Jc 


du 


t-B/2 


Hb(u) 


,  l+a'/2 

-i 
Jc 


(55) 


du  >  . 


Furthermore,  assuming  a  rectangular  window  and  applying  the  mean  value  theorem 
for  integration  produces  an  even  tighter  bound  given  by 


< 


-  sgn(o;/  +  2) 


hB(  0) 


Hb(u )  ■  mean 


'  —  B/2 


f-B/2 


Hb(u )  ■  mean 


u 

+  Jc 


U 

+  Jc 

l+a'/2 


l+a'/2 


du 


—  sgn(a/  +  2)  <  mean 

1  —B/2<u<0 


mean 

0<u<B/2 


,  u 
+  Tc 


,  U 

+  7c 

l+a'/2 


l+a'/2 


—  1 


du  >  , 


du 


(56) 
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where  the  mean  values  must  be  determined  separately  for  the  case  of  a!  =  —4,  because 
in  this  case,  the  exponent  is  -1. 

For  the  case  of  a'  ^  —4: 


mean 

—  B/2<u<0 


du 


2 

B 


'  2  fc 
4  +  a! 


2+a'/2 


—  U 


0 

-S/2 


1  [  2fc  _  2 fc  (l  _  B_\2+a'/2  _  B 

B  4  + a'  4  +  a'V  2  fcJ  2 

(57) 


mean 

0<u<B/2 


du 


2  T  2/c  /  2+0,72  B  2/e 

B  4  +  a'V  2/J  2  4  + a' 

(58) 

For  the  case  of  a'  =  —4: 


mean 

— S/2<u<0 


1  [fe 
B  J fc—B/2 


df 


n  /c 

/cj  fc—B/2 


(59) 


=  ^  [l„(/J  -  In  (/„  -  E ) 
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and 


mean 

Yi+7)"1-/ 

2 

rSc+B/ 2 

'  fc  1 

0<u<B/2 

V  fcj 

B  . 

Jfc 

J 

df 


Vc 

B 


In 


(60) 


From  Equations  (57)  through  (60)  and  recalling  from  Equation  (45)  that  the  error 
is  zero  for  the  case  of  ol  =  {—2,0},  the  bound  on  the  narrow-band  approximation 
relative  error  is 


0,  a'  €{-2,0} 


e(rect,  (3)  <  < 


!ln 


2  +  0 
2-/5 


2, 


a'  =  -4 


5 


-4sgn(a'  +  2)  [  /2  +  /3\  2+a'/2 
/5(4  +  a')  V  2  ) 


otherwise 


[  +2sgn(a/  +  2), 

(61) 

where  the  window  in  frequency  is  assumed  to  be  a  rect  function  given  by  Hpjf)  = 
rect[— B/2,  B/2],  The  bounds  for  a'  G  {  —  1,1,2}  are  shown  in  Figure  22  for  a  rect 
window.  Note  that  the  actual  error  is  well  below  this  bound  for  the  case  of  a  tapered 
window.  For  example,  the  results  from  a  numerical  analysis  of  Equation  (52)  for  a 
Hanning  window  are  shown  in  Figure  23.  In  this  case,  the  error  is  less  than  one- 
percent  for  a  fractional  bandwidth  of  0.5  or  less. 
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Figure  22.  Narrow-band  approximation  relative  error  bounds  in  Equation  (61)  as  a 
function  of  fractional  bandwidth  for  a  rectangular  window. 

5.1.3  Polar  Reformatting  Considerations. 

The  polar  format  algorithm  (PFA)  is  a  very  common  SAR  imaging  algorithm 
which  deserves  specific  consideration.  Recall  that  the  derivation  of  Equation  (47) 
assumed  an  imaging  operator  in  polar  coordinates  with  uniform  sampling  and  weight¬ 
ing  in  each  subdomain.  However,  the  PFA  features  a  rectangular  imaging  operator 
through  a  coordinate  transformation  that  eliminates  the  |/|  factor,  often  referred  to 
informally  as  the  Jacobian  [40].  Under  the  far  field  assumption  and  small  angle  ap¬ 
proximation,  0  <  20°,  described  earlier  in  Section  4.4,  the  multi-peak  model  is  still 
valid.  The  far  field  assumption  produces  a  polar  to  rectangular  transformation  given 
by  \f\dfdO  =  dfxdfy ,  /  =  sjfl  +  /2,  and  6  =  tan ~l{fy/ }x)  [76], 

Note  also  that  the  effective  sampling  density  changes  during  polar  to  rectangular 
reformatting  as  shown  in  Figure  24.  In  the  case  of  an  inscribed  rectangle  and 
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Figure  23.  Narrow-band  approximation  relative  error  in  Equation  (52)  numerically 
evaluated  as  a  function  of  fractional  bandwidth  for  a  Hanning  window. 


HxWt)  Heiff) 


(a)  Inscribed  rectangle  with  rectangular  (b)  Circumscribed  rectangle  with  polar  windows 

windows 


Figure  24.  Polar  reformatting  examples.  Circles  are  the  original  samples  in  polar 
coordinates,  and  squares  are  the  interpolated  samples  in  rectangular  coordinates.  The 
black  squares  are  assigned  a  value  of  zero. 
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rectangular  windows,  the  sample  density  and  weighting  are  uniform  as  shown  in 
Figure  24(a).  Here,  the  circles  represent  the  original  samples  in  polar  coordinates 
and  the  squares  are  interpolated  to  be  uniformly  sampled  in  rectangular  coordinates. 
As  a  result,  the  removal  of  the  Jacobian  causes  the  order  of  the  exponent  to  reduce 
from  a'  +  2  to  a'  for  Equations  (43)  through  (61). 

Alternately,  for  the  case  of  polar  windows  and  a  circumscribed  rectangle,  the 
weighting  and  zero-padding  around  the  original  annulus  varies  with  frequency,  as 
shown  in  Figure  24(b).  Here,  the  black  squares  are  assigned  a  value  of  zero.  As  a 
result,  a  Jacobian-like  weighting  of  the  phase  history  persists  even  after  the  Jacobian 
has  been  removed.  Therefore,  the  order  of  the  exponent  from  earlier  discussions  will 
remain  unchanged  due  to  the  polar  weighting  of  the  samples. 

For  hybrids  of  the  two  schemes  above,  the  exponent  will  need  to  be  adjusted 
according  to  the  effective  weighting  caused  by  the  resulting  sample  densities.  The 
polar  reformatting  algorithm  is  discussed  at  length  in  References  [24]  and  [76]. 

5. 1.3.1  Iterative  Curve  Fitting  Algorithm. 

For  the  case  of  two  subimages,  the  frequency  parameter  is  easily  obtained  from 
the  measured  subimage  pixel  intensities  using  Equation  (47).  However,  for  the  case 
of  more  than  two  subimages,  there  are  more  than  two  subimage  pixel  intensity  mea¬ 
surements.  In  this  case,  frequency  parameter  estimation  requires  curve  fitting  of  the 
measurements  to  the  set  of  normalized,  ideal  curves  expected  for  a  given  a'.  Typically, 
no  truth  data  exists  for  these  measurements;  so  any  attempt  to  normalize  the  mea¬ 
surements  creates  an  overdetermined  system  of  equations.  The  iterative  curve  fitting 
algorithm  developed  in  this  section  solves  the  overdetermined  system  of  equations 
using  a  gradient  decent  technique  [120,  10].  It  assumes  that  there  are  /  subimages 
for  the  jth  subaperture. 
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Using  the  preselected  subband  center  frequencies,  a  normalized  frequency  vector 
is  defined  as 

,,  [(/cl)“'+2,(/o2)",+2,...,(/rf)"'+2]T 

(/C)“'+2 

where  T  indicates  the  transpose  and  fc  is  the  center  frequency  of  the  full  bandwidth, 
B.  This  vector  can  be  interpreted  as  defining  the  normalized  values  (or  curves)  that  a 
set  of  ideal,  normalized  intensity  measurements  will  assume  for  a  given  a'.  Similarly, 
an  observation  vector  of  the  pixel  intensities  in  each  subimage  is  defined  as 


=  }9ij(mq,nq)\27  \g2j{mq,nq)\2, ...,  \gij(mq,  nq)\2 


21  t 


(63) 


Typically,  no  truth  data  exists,  and  so  the  correct  normalization  factor,  u,  is  unknown 
a  priori.  Therefore,  there  exists  an  overdetermined  linear  system  defined  as  ajv  = 
f(cd).  In  response,  the  iterative  curve  fitting  algorithm  obtains  an  estimate  of  ol  using 
a  greedy  search  method  that  continually  reduces  the  total  least  squares  error  of  the 
optimal  normalization  factor.  This  can  be  interpreted  as  adjusting  ol  to  find  the  best 
possible  fit  to  the  family  of  curves  defined  by  f  (cd). 

Beginning  with  Equation  (47),  the  algorithm  assigns  the  initial  value 

,  _  V \9lj\mgi  ng)\  J  _  9 

ay  /  f  \  “ 

iog  (ffl 

where  the  subscript  on  a'  indicates  the  current  iteration  of  the  algorithm.  In  each 
iteration,  the  optimal  normalization  factor  will  minimize  the  norm  of  the  residual 
expressed  as  \\cr/v  —  f(a',)||2.  This  is  determined  by  setting  the  derivative  of  the 
square  of  the  norm  of  the  residual  to  zero  [10,  120].  Thus,  the  optimal  normalization 
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factor  for  the  kth  iteration  is 


(cr1  cr) 


"k  = 


j  \  ’ 


(65) 


°'Tf  (a'k) 

where  the  '  symbol  indicates  that  this  is  an  estimate  of  the  normalization  factor.  This 
estimate  may  be  refined  by  additional  iterations  of  the  curve  fitting  algorithm. 

Next,  the  frequency  parameter  is  adjusted  by  a  scaled  version  of  the  norm  of  the 
residual  expressed  as 

5k  =  (0.95)fc||cr/4  -  f(o4)||2,  (66) 


where  the  factor  (0.95)fc  has  a  dampening  effect  that  ensures  convergence  of  the 
algorithm.  This  factor  was  selected  experimentally;  although  any  value  less  than 
unity  will  suffice,  where  higher  values  can  cause  a  slower  convergence  rate  and  lower 
values  can  cause  convergence  to  a  local  minima.  The  next  value  for  a'  is  chosen  via  a 
greedy  search  method  to  obtain  the  smallest  norm  of  the  residual.  This  is  expressed 
as 


a'k  +  5k,  \\cr/vk  ~  f(a4  +  4)||2  <  \\<r/h  -  i{a'k  -  5fc) 1 12 
°4+i  =  {  (67) 

a'k  —  Sk,  otherwise 

The  algorithm  iteratively  adjusts  a'  until  reaching  a  prescribed  amount  of  precision 
expressed  as  5k  <  0.01,  where  K  is  the  total  number  of  iterations.  At  this  point, 
the  estimate  of  the  frequency  parameter  is  finalized  as  a'  =  ot'K+x.  All  experimental 
results  contained  in  this  dissertation  used  the  exact  procedure  and  prescribed  values 
described  above. 

An  example  of  the  curve  fitting  algorithm  for  the  stable  peak  at  pixel  (1,3)  in 
Equation  (42)  of  Section  5.1.1  is  as  follows.  The  observation  vector  is 


cr  =  [17, 18,  20]r, 


(68) 


and  given  subband  center  frequencies  of  9.25  GHz,  9.5  GHz,  and  9.75  Ghz,  the  nor- 
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malized  frequency  vector  is 


f(a')  = 


[9.25Q,+2,9.5Q'+2,9.75a'+2': 

9.5Q,+2 


The  initial  value  for  a'  is  given  by 


/  =  lQg  (j) 

1  log  (If) 


-2  =  1.0871, 


and  the  initial  optimal  normalization  factor  is  given  by 


'  =  18.3095, 

<T 1  I  (Oil) 


The  adjustment  to  a1  for  this  iteration  is 


<$i  =  (0.95)1! |cr/z>i  -  fK)||2  =  0.0195. 


Because  hi  >  0.01,  another  iteration  is  required.  The  adjustment  to  a!  is 


a2  =  a\  +  hi  =  1.1066, 


because 


|cr/z>i  -f(a'i  +  hi)||2  =  0.0205  <  |  |cr/z>i  -  ^  hx)| |2  =  0.0206  (74) 
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The  following  concludes  the  curve  fitting  process: 


^  ( cr 7  <7 ) 

z>2  =  Tf(  ’  =  18.3086, 

^Tf{a2) 

S2  =  (0.95)2||<t /uq  -  fK)||2  =  0.0185. 
a3  =  a'2  +  S2  =  1.1251 
a4  =  a'3  +  S3  =  1.1426 
a'5  =  a'4  -  S4  =  1.1260 
ttg  =  olh  +  $5  =  1.1418 
a'7  =  a'6-5e  =  1.1267 


cd16  =  a'15  +  (Jig  =  1.1292  +  0.0095  =  1.1387 
Sl5  <  0.01  — *  a' =  a'16  =  1.14 

An  alternate  example  for  the  stable  peak  at  pixel  (3,2)  in  Equation  (42)  of  Section 
5.1.1  is  as  follows.  The  observation  vector  is 


<t  —  [10,20,  19]t, 


(76) 


and  the  initial  value  for  a'  is  given  by 


/ 

a1 


h*  (j) 
log  tiiH) 


10.1924. 


(77) 


Note  that  a'  =  10.1924  is  far  outside  the  expected  range  for  ideal  frequency  param¬ 
eters.  Using  Table  5  of  Section  3.2  as  a  guide,  the  ideal  range  is  a'  G  [—2,2],  As  a 
result,  the  SPLIT  algorithm  rejects  this  pixel  as  non-canonical  by  halting  the  curve 
fitting  algorithm  and  reassigning  it  a  value  of  ‘O’  in  Equation  (42)  of  Section  5.1.1. 

In  summary,  the  SPLIT  algorithm  requires  that  a  subimage  peak  meet  two  criteria 
before  extracting  the  frequency  parameter.  First,  the  peak  must  be  stable,  and  sec- 
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ond,  it  must  have  a  frequency  parameter  that  is  within  or  nearby  the  expected  range 
for  canonical  scatterers.  For  the  experimental  results  contained  in  this  dissertation 
the  estimated  frequency  parameter  was  required  to  be  within  the  prescribed  ranges 
of  a\  G  [—6,6]  and  a'K+l  G  [—4,4], 

5.1.4  Polarimetric  Parameter  Estimation. 

The  Krogager  scattering  matrix  parameters  are  estimated  from  fully-polarized 
SAR  data.  This  requires  the  availability  of  at  least  three  phase  histories,  two  co¬ 
polarization  phase  histories  (HH  and  VV)  and  one  cross-polarization  phase  history 
(HV  or  VH).  The  Krogager  parameters  were  introduced  earlier  and  discussed  at  length 
in  Sections  2.2.3  and  2.3.3. 

The  Krogager  parameters,  and  ne,  for  a  subimage  pixel  having  a  stable  peak 
are  easily  obtained  by  a  coherent  summation  of  the  pixel  values  from  the  polarization- 
diverse  subimages.  The  parameters  represent  the  proportion  of  odd-bounce  and  even- 
bounce  scattering  in  each  pixel  of  the  subimage,  where  it  is  convenient  to  define  the 
proportion  of  helical  scattering  as  Kh  =  1  -  -  ne. 

In  general,  it  has  been  noted  that  the  polarimetric  parameters  of  canonical  scat¬ 
terers  can  be  estimated  with  greater  accuracy  than  the  other  parameters.  For  this 
reason,  the  latest  variant  of  the  canonical  scatterer  classifier  developed  at  The  Ohio 
State  University  estimates  canonical  scatterer  parameters  from  an  image  segment  us¬ 
ing  a  tiered  approach,  rather  than  a  joint  estimation  approach.  Because  estimates  of 
the  polarization  parameters  are  assumed  to  have  the  highest  statistical  confidence, 
these  are  estimated  first,  so  as  to  restrict  the  possible  values  of  the  other  parameters 
in  subsequent  tiers  of  the  estimation  scheme  [73]. 

Deference  to  the  polarimetric  parameters  is  also  adopted  in  this  dissertation.  The 
assumption  is  that  the  frequency  parameter,  a',  has  a  much  higher  variance  than 
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the  Krogager  parameters  and  is  much  more  sensitive  to  modeling  error,  signal  noise, 
and  clutter.  Because  of  this  assumption,  some  ambiguities  in  the  feature  space  are 
resolved  by  favoring  instances  of  strong  polarimetric  response.  For  example,  according 
to  Table  5  from  Section  3.2,  the  ideal  feature  vector  for  a  horizontal  cylinder  is 
wa  =  [a',K0,Ke]  =  [—1,1,0],  while  the  ideal  feature  vector  for  a  horizontal  wire  is 
w6  =  [-2,  0.5, 0.5], 

Consider  the  case  of  an  extracted  feature  vector,  v  =  [—2,1,0],  where  the  Eu¬ 
clidean  norms  are  |jwa  —  v 1 1 2  =  1  and  ||w&  —  v 1 1 2  =  0.707.  Clearly  the  norm  is  smaller 
for  the  horizontal  wire;  however,  ‘horizontal  cylinder’  is  the  best  classification  deci¬ 
sion  because  both  norms  are  large  and  the  polarimetric  response  better  matches  that 
of  the  horizontal  cylinder.  In  this  way,  the  concept  of  polarimetric  preference  will 
impact  the  partitioning  of  the  feature  space  for  the  SPLIT  algorithm,  as  described 
later  in  Section  5.3. 

5.2  Feature  Vector  Summation 

The  SPLIT  algorithm  extracts  the  frequency  parameter  from  each  co-polarization 
channel.  The  result  is  a  single  pixel  having  two  estimates  of  the  frequency  parame¬ 
ter,  a1,  when  both  HH  and  VV  polarizations  are  available.  In  addition,  the  SPLIT 
algorithm  extracts  the  polarization  parameters  from  each  subimage.  This  results  in 
a  single  pixel  having  /  estimates  of  the  polarimetric  parameters,  [k0,  Ke].  Finally,  the 
SPLIT  algorithm  extracts  frequency  and  polarimetric  features  for  each  subaperture. 
Therefore,  it  is  possible  for  the  same  pixel  to  be  assigned  multiple  feature  vectors, 
each  extracted  from  a  different  subaperture.  In  these  instances,  the  feature  vectors 
for  a  single  pixel  are  combined  using  a  weighted  average.  Recall  that  pixels  failing 
the  canonical  scatterer  criteria  of  the  SPLIT  algorithm  are  assigned  a  feature  vector 
weight  of  zero. 
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For  the  case  of  both  HH  and  VV  polarizations,  the  frequency  parameter  is  equal 


to  the  weighted  average 


ajpajp 


ajp 


(78) 


p 

where  <Jjp  =  minj{| gijP(mq,  nq)  |2}  is  the  minimum  subimage  pixel  intensity  for  a  given 
subaperture  and  polarization. 

Likewise,  the  polarimetric  parameters  are  equal  the  to  weighted  vector  sum 


^  '  &ij  K0,  K-e  ij 


[l^o  5  ^e\j 


E 


a, 


(79) 


where  atJ  =  minp{| gijP(mq,  nq)  |2}  is  the  minimum  subimage  pixel  intensity  for  a  given 
subband  and  subaperture. 

Finally,  referring  to  Equation  (41),  the  feature  vector  for  a  single  pixel  observed 
at  multiple  subapertures  is  equal  to  the  weighted  vector  sum 


v(mq,nq )  = 

j 


(80) 


where  <jj  =  nq)\2}  is  the  minimum  subimage  pixel  intensity  for  a  given 

subaperture.  This  vector  is  referred  to  as  the  average  feature  vector. 

An  example  of  cross-aperture  feature  vector  summation  is  as  follows.  Consider 


the  4x4  feature  vector  arrays  for  subapertures  j  —  1  and  j  =  2  given  as 


V1  = 


0  0  0  0 

0  0  0  0 

17[1.14,  0.22,  0.67]  0  0  0 

0  0  0  0 


V,  = 


0 

0 

14[0.99,  0.25,0. 

0 


0  42[1.79, 0.82, 0.11]  0 

0  0  0 

0  0  0 

0  0  0 


(81) 


(82) 


where  0  is  the  zero  vector.  Note  that  the  feature  vector  array  in  Equation  (81)  is 
derived  from  the  earlier  examples  in  Equations  (42)  and  (75).  The  weighted  vector 
sum  from  Equation  (80)  computed  for  each  pixel  results  in  an  array  of  average  feature 
vectors  expressed  as 


0  0  [1.79,0.82,0.11]  0 

0  0  0  0 

[1.07,0.23,0.63]  0  0  0 

0  0  0  0 


(83) 


where  the  pixels  with  non-canonical  scattering  have  weights  of  zero  and  contribute 
nothing  to  the  vector  sum. 


5.3  Canonical  Scatterer  Classification 

The  feature  space  for  canonical  scatterers  can  be  partitioned  using  the  ideal  fea¬ 
ture  vectors  in  Table  5.  However,  these  do  not  uniformly  fill  an  entire  Euclidean 
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feature  space.  Therefore,  additional  ideal  feature  vectors  are  added  with  deference 
to  the  polarimetric  parameters,  as  described  earlier  in  Section  5.1.4.  First,  the  heli¬ 
cal  scatterer  class  is  assigned  five  ideal  feature  vectors,  [—2,0,0],  [—1,0,0],  [0,0,0], 
[1,  0,  0],  and  [2,  0,  0].  The  assumption  is  that  helical  scattering  can  be  associated  with 
any  value  of  a'.  Next,  the  dihedralo  scatterer  class  is  assigned  three  ideal  feature  vec¬ 
tors,  [0,0,1],  [—1,0,1],  and  [—2,0,1].  Last,  the  cylinder0  scatterer  class  is  assigned 
two  ideal  feature  vectors,  [—1, 1,  0]  and  [—2, 1,  0]. 

The  result  is  a  set  of  seventeen  ideal  feature  vectors,  IV,  assigned  within  a  set 
of  ten  scatterer  classes,  C.  Scatterer  classification  is  accomplished,  pixcl-by-pixel,  by 
comparing  each  non-zero  feature  vector  in  the  array  of  average  feature  vectors,  V,  to 
the  set  of  ideal  feature  vectors,  W.  The  classification  decision  can  be  accomplished 
using  a  likelihood  ratio  test  informed  by  the  statistical  properties  of  each  class  and 
parameter,  if  available.  However,  for  illustrative  purposes,  a  simple  a  least  squares 
classifier  is  preferred  for  demonstrating  the  SPLIT  algorithm.  A  least  squares  classifier 
assumes  uniform  costs  for  misclassification  and  equal  prior  probabilities  of  each  class. 
The  least  squares  classifier  and  an  associated,  non-statistical  measure  of  fitness  is 
presented  in  this  section  and  used  throughout  the  remainder  of  this  dissertation. 

5.3.1  Least  Squares  Classifier. 

For  least  squares  classification,  each  non-zero  average  feature  vector  is  compared 
to  the  set  of  seventeen  ideal  feature  vectors,  W,  using  the  Euclidean  norm.  Hence, 
each  subimage  pixel  is  declared  a  member  of  the  class  associated  with  the  ideal  feature 
vector  for  which  the  Euclidean  norm  is  a  minimum.  The  feature  vector  decision  is 
succinctly  represented  as 


d  =  argrnin  \\Wk  -  v||2, 


(84) 


A 

trihedral 

0 

dihedralgo 

□ 

cylindergQ 

V 

top  hat 

O 

sphere,  plate 

* 

edge/wire9Q 
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dihedral0 

■ 

cylinder0 

* 

edge/wireg 

<1 

helical 

Figure  25.  Feature  space  partitioned  into  scatterer  classes  according  to  a  least  squares 
classifier. 


where  k  G  {1,2, ...,  17}  provides  an  ordinal  list  of  the  ideal  feature  vectors  and  d  is 
the  index  of  the  feature  vector  leading  to  the  smallest  norm  in  Equation  (84).  Hence, 
Equation  (84)  divides  the  feature  space  into  seventeen  regions  according  to  a  least 
squares  “best  fit”  for  the  ideal  canonical  scatterers,  as  shown  in  Figure  25. 

Finally,  the  scatterer  classification  decision  is  expressed  using  a  scatterer  classifi¬ 
cation  operator,  S,  as 

<?{v}=C(Wd),  (85) 

where  the  class  associated  with  ideal  feature  vector  determines  the  assigned  class. 
For  example,  least  squares  classification  using  the  array  of  average  feature  vectors  in 


Equation  (83)  results  in 


—  trihedral  — 


V 


tophat  — 


(86) 


This  is  due  to  the  fact  that  the  average  feature  vectors  v(3, 1)  =  [1.79,  0.82,  0.11]  and 
v(l,3)  =  [1.07,0.23,0.63]  are  closest  to  the  ideal  feature  vectors  [2,1,0]  and  [1,0,1], 
which  in  turn,  are  associated  with  the  trihedral  class  and  top  hat  class,  respectively. 


5.3.2  Measure  of  Fitness. 


Some  average  feature  vectors  may  reside  near  a  boundary  between  two  regions. 
Therefore,  it  is  instructive  to  provide  an  additional  measure  of  fitness.  The  measure 
of  fitness  is  defined  as 


7=1- 


II Wd  ~  v||2 


min  1 1 W4  —  v  1 1 2  ’ 
yvkic(wd) 


(87) 


where  the  denominator  is  the  Euclidean  norm  for  the  next  best  fit  for  ideal  feature 
vectors  not  associated  with  the  declared  class,  C  (Wb).  Using  a  ratio  of  Euclidean 
norms  causes  the  measure  of  fitness  to  be  in  the  range  7  G  [0, 1]. 

For  example,  beginning  with  Equation  (83),  the  measures  of  fitness  are 


||[1, 0, 1]  -  [1.07, 0.23, 0.63]||2 
II [1,0,0]  -  [1.07,  0.23,  0.63]|| 


(88) 


and 


||[2,1,0]~  [1.79, 0.82, 0.11]||2 
||[1, 1,0]  -  [1.79,  0.82,  0.11]||2 


0.69, 


(89) 


where  the  next  best  fits  are  a  helical,  [1,  0,  0],  and  cylinder,  [1, 1,  0],  respectively.  Thus, 
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the  array  of  measures  of  fitness  for  the  example  in  Equation  (83)  is 


-  -  0.69  - 


r 


0.35  -  -  - 


(90) 


where  Equations  (86)  and  (90)  represent  example  classification  results  produced  by 
the  SPLIT  algorithm. 


5.4  Experimental  Results 

This  section  presents  experimental  results  using  measured  and  simulated  SAR 
data.  The  complexity  of  the  scene  increases  with  each  experiment.  First,  isolated 
canonical  scatterers  are  classified  to  demonstrate  and  validate  the  usefulness  of  the 
multi-peak  model  for  scatterer  classification.  Then,  inverse  SAR  (ISAR)  radar  mea¬ 
surements  of  a  smooth  metallic  body,  called  D2,  are  evaluated.  Even  though  the 
measurements  are  not  fully-polarimetric,  the  D2  provides  an  excellent  case  study  for 
modeling  a  simple  object  as  a  compilation  of  canonical  scatterers.  Next,  simulated 
data  of  civilian  vehicles  provide  a  more  challenging  scenario  for  scatterer  classification. 
These  illustrate  both  the  benefits  and  limitations  of  scatterer  classification  by  phase 
domain  decomposition  for  complex  objects.  Finally,  the  Gotcha  multipass  data,  fea¬ 
turing  measured  SAR  data  of  a  calibration  targets,  is  examined.  Because  the  Gotcha 
data  has  only  moderate  bandwidth,  poor  coherency  between  pulses,  and  high  clutter 
energy,  the  classification  results  are  marginal,  as  expected. 

Arrays  similar  to  those  in  Equations  (86)  and  (90)  are  used  to  create  an  overlay 
to  display  the  classification  results  onto  the  SAR  image.  The  classification  results 
are  displayed  using  the  symbols  and  colors  for  the  ideal  feature  vectors  in  Figure  25, 
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while  the  measure  of  fitness  is  displayed  by  symbol  size.  Note  that  the  helical  scatterer 
is  given  a  neutral  color  in  the  legend,  but  each  of  the  five  ideal  feature  vectors  will 
actually  give  a  unique  color  in  the  overlay  for  each  value  of  a! .  The  colors  in  ascending 
order  are  magenta  (cd  =  —2),  red  (a!  =  —1),  yellow  ( a '  =  0),  green  ( a '  =  1),  and 
cyan  ( a '  =  2). 

The  images  presented  in  this  Section  were  produced  by  integrating  the  SPLIT- 
based  scatterer  classification  algorithm  with  a  fast  convolution  backprojection  imaging 
algorithm  as  described  in  the  next  chapter.  The  integrated  algorithm  employs  domain 
decomposition  to  produce  coarse-resolution  subimages  for  subsequent  analysis  by  the 
scatterer  classifier.  Then  it  interpolates,  weights,  and  sums  the  subimages  to  form  a 
fine-resolution  SAR  image. 

5.4.1  Accuracy  of  the  Multi-peak  Model. 

It  is  desirable  to  succinctly  present  the  accuracy  of  the  new  multi-peak  model  for 
all  canonical  scatterer  types  and  at  several  non-zero  orientation  angles.  In  addition,  it 
is  helpful  to  illustrate  how  the  multi-peak  model  in  Equation  (37)  is  useful  for  classi¬ 
fication  of  canonical  scatterers  in  SAR  imagery.  Consequently,  this  section  illustrates 
the  applicability  of  the  model  to  scatterer  classification  while  reiterating  some  of  the 
benefits  and  limitations  of  the  model  which  were  discussed  previously. 

Because  frequency  parameter  estimation  in  Equation  (47)  uses  the  new  multi¬ 
peak  model  parameterized  by  a',  it  is  instructive  to  re-evaluate  the  applicability  of 
the  multi-peak  model  in  Equation  (37)  for  use  with  the  SPLIT  algorithm.  From 
Equation  (31)  and  assuming  9C  =  0  without  loss  of  generality,  the  ratio  of  integrated, 
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Ratio  of  Integrated,  Windowed  Sine  Functions 


0  (degrees) 


1 

0.98 

0.96 

0.94 
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0.9 


Figure  26.  Under  wide-angle  conditions,  the  ratio  of  integrated,  windowed  sine  func¬ 
tions  equals  the  inverse  ratio  of  sine  frequencies.  Here,  fi  =  10-Ghz,  /2  =  9-GHz,  and 

e0  =  0°. 


windowed  sine  functions  must  satisfy  the  relationship 


j  Hq(9)  sine  (f/iLsin(6>  -  0O))  j 

He(9)  sine  f2Lsm(9-90))d9  h 


(91) 


where  0  >  Qmin  and  the  minimum  aperture  width,  @m,n,  can  be  determined  by 
numerical  integration.  Numerical  results  for  Equation  (91)  using  a  Hanning  window 
are  shown  in  Figure  26,  where  f\  =  10  GHz  and  /2  =  9  GHz.  The  graph  on  the  right 
indicates  that  for  0  >  10°  the  relationship  in  Equation  (91)  is  satisfied,  where  f2/fi  = 
0.9.  Earlier  analyses  indicate  that  this  relationship  holds  for  varying  frequencies  and 
orientations  angles  as  long  as  the  aperture  width  is  restricted  to  0  >  10°  and  the 
window  function  in  azimuth  is  tapered. 

Simulated  phase  histories  were  produced  for  fifteen  canonical  scatterers  listed  in 
Table  8  using  canonical  scattering  models  from  Reference  [75].  This  data  set  contains 
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Table  8.  Simulated  Canonical  Scatterer  Specifications. 


Index 

Description 

Height 

(m) 

Width 

(m) 

Radius 

(m) 

L 

(m) 

T 

(deg) 

(A) 

square  trihedral 

1 

1 

n/a 

0 

0 

(B) 

dihedralgo 

1 

1 

n/a 

0 

90 

(C) 

cylinder^ 

1 

n/a 

0.5 

0 

90 

(D) 

top  hat 

1 

n/a 

0.5 

0 

0 

(E) 

point  scatterer 

n/a 

n/a 

n/a 

0 

n/a 

(F) 

sphere 

n/a 

n/a 

0.5 

0 

n/a 

(G) 

plate,  15A 

1 

0.5 

n/a 

0.5 

0 

(H) 

plate,  30A 

1 

1 

n/a 

1 

0 

(I) 

plate,  60A 

1 

2 

n/a 

2 

0 

(J) 

dihedralo,  15A 

1 

0.5 

n/a 

0.5 

0 

(K) 

dihedralo,  30A 

1 

1 

n/a 

1 

0 

(L) 

dihedralo,  60 A 

1 

2 

n/a 

2 

0 

(M) 

cylinder0, 15A 

n/a 

0.5 

0.5 

0.5 

0 

(N) 

cylinder0,  30 A 

n/a 

1 

0.5 

1 

0 

(0) 

cylinder0,  60 A 

n/a 

2 

0.5 

2 

0 

Table  9.  Radar  Measurement  Data  Parameters. 


Parameter 

Value 

Polarizations 

HH,  VV,  VH 

LFM  Frequencies 

8.6061  :  0.0081  :  10.5939  GHz 

Azimuth 

-5.0000  :  0.0859  :  4.9695  degrees 

Elevation  Angle 

45  degrees 

specific  examples  of  the  ideal  canonical  scatters  in  Table  5,  to  include  six  different 
point  scatterers  and  three  different  distributed  scatterers  of  varying  lengths.  The 
lengths  were  specifically  chosen  to  illustrate  the  limitations  of  the  high  frequency  ap¬ 
proximation  at  X-band.  The  phase  histories  were  simulated  using  a  circular  aperture 
according  to  the  specifications  in  Table  9.  Here,  the  azimuth  angle  is  referenced 
counter-clockwise  from  the  negative  y-axis,  and  the  frequency  and  azimuth  are  listed 
as  start :  increment :  end  values. 

Classification  accuracy  is  dependent  upon  model  accuracy,  bandwidth,  and  inter¬ 
ference  from  neighboring  canonical  scatterers,  clutter,  and  noise.  However,  in  order 
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Figure  27.  Classification  results  for  0  =  10°.  Scatterers  (A)  to  (F)  are  point  scatterers 
and  (G)  to  (O)  are  distributed  scatterers  of  varying  lengths. 


to  isolate  the  effect  of  model  accuracy  on  classification  accuracy,  the  simulations  are 
free  of  interference.  The  images  and  classification  results  below  were  produced  using 
the  integrated  algorithm  with  1  =  3  subbands.  Also,  the  images  are  truncated  to 
show  only  the  top  45  dB  of  image  intensity. 

For  the  case  of  0  =  10°  shown  in  Figure  27,  a  single  aperture  is  used  to  produce 
the  classification  results  and  images.  Here,  the  classification  results  are  annotated  on 
the  image  using  the  symbology  from  Figure  25.  For  this  case,  all  classification  results 
are  correct  because  the  multi-peak  model  is  accurate  when  the  wide-angle  condition 
is  satisfied  and  all  distributed  scatterers  are  greater  than  ten  wavelengths  in  length. 
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Figure  28.  Classification  results  for  narrow-angles.  For  the  distributed  scatterers, 
accuracy  decreases  as  0  or  L  decreases. 


In  this  case,  the  longer  distributed  scatterers  (I),  (L)  and  (0),  each  produce  a  group 
of  in-line  peaks.  The  in-line  peaks  are  a  manifestation  of  Gibbs  phenomenon  due 
to  the  imaging  operator,  and  illustrate  a  characteristic  of  distributed  scatterers  from 
which  the  multi-peak  model  derives  its  name. 

For  the  case  of  0  =  5°  shown  in  Figure  28(a),  three  evenly  spaced,  overlapping 
subapertures  were  used  to  produce  the  classification  results  and  images.  In  this  case, 
the  classification  results  are  correct  for  all  point  scatterers  (A)  to  (F)  and  for  the 
longer  distributed  scatterers,  as  expected.  Also  as  expected,  the  classification  results 
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Figure  29.  Numerical  results  for  Eq.  (31)  where  Hq  is  a  Hanning  function  with  region 
of  support  equal  to  0.  The  distributed  scatterer  is  of  length  L  =  0.5-meter. 

are  incorrect  for  the  shortest  distributed  scatterers,  (G),  (J),  and  (M),  because  the 
accuracy  of  the  new  multi-peak  model  is  degraded  for  the  case  of  0  =  5°  and  L  =  15A. 
These  results  were  predicted  by  numerical  analysis,  where  Figure  26  reveals  that  the 
model  is  only  accurate  when  L  <  25 A  for  the  case  of  0  =  5°.  The  two  peaks  which 
were  correctly  classified  came  from  analysis  of  the  subapertures  at  9q  =  [—2.5,  2.5]. 
It  can  be  seen  in  Figure  29  that  the  model  is  valid  at  0  =  5°  because  A  is  essentially 
independent  of  frequency  for  this  case.  However,  the  value  at  A  is  very  small  when 
0  =  5°,  which  results  in  low  peak  intensities.  Low  peak  intensities  are  a  well-known 
characteristic  of  long  scatterers  imaged  at  off-broadside  aspects  [83] . 

For  the  case  of  0  =  2.5°  shown  in  Figure  28(b),  seven  evenly  spaced,  overlapping 
subapertures  were  used  to  produce  the  classification  results  and  images.  In  this 
case,  the  classification  results  are  correct  for  all  point  scatterers  (A)  to  (F)  and  for 
the  longest  distributed  scatterers,  as  expected.  Also  as  expected,  the  classification 
results  are  incorrect  for  the  shorter  distributed  scatterers,  (G),  (H),  (J),  (K),  (M), 
and  (N),  because  the  accuracy  of  the  new  multi-peak  model  has  further  degraded  for 
this  case.  The  reasons  for  this  are  analogous  to  the  case  of  0  =  5°,  and  the  results 
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Table  10.  Peak  Intensity  Classification  Results,  0  =  10°. 


Index 

(  x 

,  y  ) 

(  x 

>  y  ) 

(  y  ) 

(  x 

■  y  ) 

(  y  ) 

ec 
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(7 
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a 

ex'  a 
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a 
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- 

0.0° 

- 

- 
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- 
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- 

0.0° 

- 

- 
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- 

-  1  - 

(C) 
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- 

0.0° 

- 

- 
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- 
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- 

-  1  - 
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- 

0.0° 

- 

- 

0.00  I  0  dB 

- 

-  1  - 
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- 

0.0° 

- 

- 
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- 

-  1  - 
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(0,0) 

- 
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- 

- 

0.01  I  0  dB 

- 
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- 
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- 

- 
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- 
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- 
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- 
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(0,-0. 76) 

(0.76,-0.76) 

0.0° 

-0.94 

0  dB 

- 

-0.99  |  0  dB 

- 

-0.94  |  0  dB 

are  as  expected  based  on  Figure  26.  Evidently,  the  model  is  accurate  for  L  =  60A 
when  0  =  2.5°,  as  can  be  seen  by  examining  (I),  (L),  and  (O). 

If  desired,  the  classification  process  for  this  experiment  can  be  verified  by  exam¬ 
ining  the  extracted  frequency  parameters  provided  in  Tables  10,  11,  and  12.  All 
intensity  values  are  normalized  to  the  highest  intensity  sample  for  each  scatterer, 
individually.  The  Krogager  parameters  are  not  listed,  because  the  simulated  data 
produced  almost  ideal  Krogager  values  in  each  case. 

In  summary,  all  classification  results,  including  incorrect  classifications,  were  nicely 
predicted  using  the  multi-peak  model.  The  classification  results  were  correct  for  the 
canonical  point  scatterers  (A)  to  (F)  and  the  longer  distributed  canonical  scatterers 


Table  11.  Peak  Intensity  Classification  Results,  0  =  5°. 
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(I),  (L),  and  (0),  in  every  case.  For  the  incorrect  classifications,  the  experiment  il¬ 
lustrates  how  model  accuracy  decreases  gradually  as  0  or  L  decreases.  These  results 
indicate  that  the  multi-peak  model  can  be  used  to  classify  distributed  scatterers  by 
comparing  the  peak  intensities  of  subimages  produced  by  phase  history  decomposi¬ 
tion.  However,  the  aperture  should  be  large  (0  >  10°),  particularly  for  distributed 
scatterers  of  shorter  physical  length. 
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(a)  Frontal  aspect.  (b)  Rear  aspect.  (c)  Dimensions. 

Figure  30.  D2  airframe  [4], 


5.4.2  D2  Measurements. 

This  section  presents  experimental  results  obtained  from  measurements  of  a  smooth 
metallic  body  called  D2.  The  measured  data  was  obtained  from  the  AFIT  indoor  RCS 
measurement  range  described  in  Appendix  A.  Photographs  and  dimensions  of  the  D2 
are  provided  in  Figure  30,  and  a  photo  of  the  measurement  configuration  is  provided 
in  Figure  31.  Additional  information  regarding  D2  and  its  measured  data  can  be 
found  in  Reference  [4]. 

The  parameters  for  an  ISAR  measurement  are  given  in  Table  13,  and  a  2D  im¬ 
age  overlaid  with  classification  results  combined  from  both  HH  and  VV  channels  is 
shown  in  Figure  32.  Note  how  the  straight  edges  of  the  D2  are  correctly  classified  as 
horizontal  edges  and  how  the  D2-pedestal-D2  interactions  are  consistently  classified 
as  trihedrals  (or  horizontal  dihedrals).  These  correct  classifications  occur  in  spite 
of  the  fact  that  these  measurements  are  in  the  near  field,  which  violates  one  of  the 
assumptions  of  the  multi-peak  model. 

The  near  field  distortion  creates  an  artificial  curvature  that  may  explain  the  ‘cylin¬ 
der’  classifications  along  the  back  edge  near  (x  =  0.9,  y  =  0.1)  m  in  Figure  32.  Like- 
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Figure  31.  Photo  of  D2  airframe  from  a  rear  aspect  [4].  The  D2  is  on  the  target 
pedestal  at  the  AFIT  Indoor  RCS  measurement  range  with  the  radar  antennas  shown 
in  the  background. 


Table  13.  Experimental  Parameters  for  ISAR  Measurements  of  D2. 


Parameter 

Value 

Polarizations 

VV  and  HH 

LFM  Frequencies  (GHz) 

6.528  :  0.023  :  18.005 

Fractional  bandwidth 

/3  =  0.936 

Elevation 

0° 

Azimuth 

0.2°  :  0.2°  :  360.0° 

Azimuth  subwindow 

10°  Hanning 

Frequency  subwindow 

halfband  Hanning 

Number  of  subbands 

1  =  5 

Intensity  threshhold 

top  25  db 

Subaperture  summation 

coherent 

Oversampling 

none 

wise,  the  curved  edges  near  (0.5,  ±0.75)  m  are  likely  classified  as  cylinders  for  similar 
reasons.  As  described  in  Section  2.2.1,  it  has  been  shown  that  curved  edges  create 
point  scatterers  having  a  frequency  response  parameterized  by  a  =  —1  [62,  84,  121]. 
The  curved  edge  canonical  scatterer  is  not  included  in  the  SPLIT  algorithm  because 
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Figure  32.  2D  ISAR  measurement  with  superimposed  HH  and  VV  images  and  classi¬ 
fication  data. 


its  low  RCS  makes  it  likely  to  be  unresolved  in  most  SAR  applications.  However,  for 
low  RCS  airframes  where  smaller  RCS  mechanisms  are  prevalent,  the  ideal  feature 
vector  [a' ,  n0,  /te]  =  [—1,  0.5,  0.5]  could  easily  be  included  to  account  for  curved  edges. 

Unfortunately,  it  is  not  entirely  clear  why  the  front  tip  at  (—1,0)  m  is  classified 
as  a  sphere  or  other  scatterer  with  a 1  =  0.  It  may  be  due  to  a  traveling  wave  down 
the  long  edge  [83].  Note  also,  that  interference  from  neighboring  scatterers  in  the 
center  of  the  pedestal  seems  to  have  caused  some  confusion  in  the  classifier.  This 
interference  is  evident  by  the  random  classification  results  within  a  circle  of  radius 
0.1  m  in  the  center  of  the  image.  Interference  from  the  pedestal  is  also  likely  to  have 
impacted  the  classification  results  near  (0.1,  0.5)  m.  Overall,  the  classification  results 
in  Figure  32  could  prove  useful  to  an  analyst,  particularly  the  trihedral  and  edge 
classifications  which  appear  symmetric.  It  is  also  likely  that  an  experienced  analyst 
will  be  able  to  disregard  classification  results  in  areas  deemed  likely  to  produce  high 
levels  of  interference  based  on  knowledge  of  the  image  context. 
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Figure  33.  Facet  model  of  Toyota  Tacoma  2-door  from  120-degrees  azimuth  and  30- 
degrees  elevation. 


Table  14.  Experimental  Parameters  for  Toyota  Tacoma  Data  Dome. 


Parameter 

Value 

Polarizations 

VV,  HH,  and  HV 

Bandwidth  (GHz) 

6.9226  :  0.01046  :  12.2774 

Fractional  bandwidth 

(3  =  0.558 

Elevation 

45° 

Azimuth 

0.0625°  :  0.0625°  :  360.0° 

Azimuth  subwindow 

10°  Hanning 

Frequency  subwindow 

halfband  Hanning 

Number  of  subbands 

1  =  5 

Intensity  threshhold 

top  40  db 

Subaperture  summation 

multilook 

Oversampling 

none 

5.4.3  Civilian  Vehicle  Data  Domes. 

The  SPLIT  classification  algorithm  was  also  run  for  a  simulated  data  set  of  civilian 
vehicles,  which  are  structurally  more  complex  than  D2.  A  facet  model  of  the  Toyota 
Tacoma  vehicle  from  the  Air  Force  Research  Laboratory  civilian  vehicle  (CV)  data 
domes  [124]  is  shown  in  Figure  33  and  the  classification  results  are  given  in  Figure 
34.  The  associated  data  parameters  are  listed  in  Table  14.  The  truck  bed, 
in  particular,  provides  a  scatterer  rich  environment,  creating  interference  to  stress 
the  classification  algorithm.  In  addition,  an  the  elevation  angle  of  45°  combined 
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Figure  34.  Annotated  SAR  image  of  a  simulated  Toyota  Tacoma  pointed  in  the  — x 
direction. 


with  ground  plane  imaging  was  chosen  so  that  layover  would  separate  some  of  the 
prime  scattering  mechanisms  from  potential  sources  of  interference.  Generally,  the 
primary  scattering  mechanisms  on  civilian  vehicles  are  cylindrical  returns  from  the 
top  edges  of  the  vehicle  and  dihedral  returns  from  the  sides  of  the  vehicle  and  the 
ground,  as  depicted  in  Figure  35  [44,  49].  Figure  36  shows  the  cylinder0  and  dihedralo 
classifications  separately  in  Figures  36(a)  and  36(b),  respectively,  to  highlight  these 
primary  scattering  mechanisms. 

Despite  the  complexity  of  the  target  and  the  simplicity  of  the  scatterer  classifier, 
three  dominant  scattering  mechanisms  are  identifiable  as  well  classified.  First,  the 
cylindrical  edges  of  the  cab  roof  produce  a  large  circular  footprint  with  a  radius  of 
about  2.3  nr  centered  at  point  (—0.3,  0)  nr.  The  large  circle  appears  in  the  image 
due  to  layover  effect  [76].  Most  of  the  image  peaks  associated  with  this  mechanism 
are  correctly  classified  as  a  cylinder  with  a  tilt  angle  of  zero  degrees.  Second,  the 
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Figure  35.  The  dominant  backscatter  mechanisms  for  passenger  vehicles  are  a  single 
bounce  from  the  top  edge  (solid  lines)  and  double  bounce  from  the  dihedral  formed 
with  the  ground  plane  (dotted  lines). 


>> 


-4  -3  -2-10  1 

x  (meters) 


-1  0  1 
x  (meters) 


(a)  Horizontal  cylinders  only  (b)  Horizontal  dihedrals  only 

Figure  36.  Classification  results  for  the  Toyota  Tacoma  for  specific  scatterers. 


right  angle  formed  between  the  rear  bumper  and  tailgate  produces  a  vertical  line 
about  1.5  nr  long  centered  at  point  (2.9,0).  The  displacement  from  the  tailgate 
located  at  x  ~  2.2  nr  is  due  to  layover  effect.  Most  of  the  inrage  peaks  associated 
with  this  mechanism  are  correctly  classified  as  a  dihedral  with  a  tilt  angle  of  zero 
degrees.  Third,  the  right  angles  formed  between  the  truck  body  and  the  ground  plane 
produce  a  rectangular  footprint  about  1.5  nr  along  the  y-axis  and  4.5  nr  along  the 
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Figure  37.  Annotated  SAR  images  of  three  automobiles,  each  pointed  in  the  —x  direc¬ 
tion. 


(a:)-axis  centered  at  (0,0).  Many  of  the  image  peaks  associated  with  this  mechanism 
are  correctly  classified  as  a  dihedral  with  a  tilt  angle  of  zero  degrees;  however,  a 
preponderance  of  clutter,  particularly  from  the  wheels  and  wheel  wells  compromise 
the  classification  results.  All  three  mechanisms  are  examples  of  distributed  canonical 
scatterers  found  in  the  geometry  of  the  target  which  are  well  classified  using  the 
multi-peak  model.  These  scattering  mechanisms  are  also  identifiable  for  the  sedans 
featured  in  Figure  37,  where  an  inner  rectangle  of  horizontal  dihedrals  is  ringed  by 
horizontal  cylinders. 
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Table  15.  Experimental  Parameters  for  Gotcha  Data  Set. 


Parameter 

Value 

Polarizations 

HH,  VV,  HV,  and  VH 

LFM  Frequencies  (GHz) 

9.2881  :  0.0014  :  9.9105 

Fractional  bandwidth 

(3  =  0.0648 

Mean  Elevation 

43° 

Ground  Range  Resolution 

0.17  meters 

Azimuth 

-179.9934°  :  0.0084°  :  179.9982° 

Azimuth  subwindow 

10°  Hanning 

Frequency  subwindow 

halfband  Hanning 

Number  of  subbands 

1  =  3 

Intensity  threshhold 

top  45  db 

Subaperture  summation 

coherent 

Oversampling 

2x 

Autofocus 

yes,  given  in  data  set 

5.4.4  Gotcha  Public  Release  Data. 

In  2007,  the  Air  Force  Research  Laboratory  publicly  released  a  challenge  data 
set  featuring  several  circular  SAR  measurements  of  a  parking  lot  scene  [25].  One  of 
these  measurements,  called  pass  number  six,  appeared  to  have  superior  stability  in 
the  circular  aperture  as  evidenced  by  images  with  better  focus.  The  radar  data  from 
pass  number  six  is  used  exclusively  in  the  following  analysis,  and  its  parameters  are 
listed  in  Table  15.  Note  that  the  LFM  Frequency  and  Azimuth  values  are  listed  as 
start :  increment :  end  values  from  the  Gotcha  data. 

Ideally,  the  SAR  system  collects  data  at  regular  intervals  along  the  circular  flight 
path  so  that  the  spectral  domain  is  uniformly  sampled  in  azimuth.  However,  this 
is  generally  not  the  case.  Flight  dynamics  and  variable  winds  cause  the  ground 
speed  and  altitude  of  the  airborne  radar  to  vary  throughout  the  flight  path.  These 
conditions  are  controlled  to  the  best  extent  possible,  resulting  in  azimuth  sampling 
that  is  pseudo-uniform  and  an  elevation  angle  which  varies  slightly  around  the  mean. 
Such  variances  are  typical  in  circular  SAR  data  and  are  not  expected  to  adversely 
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Figure  38.  Gotcha  calibration  target  types  and  locations  taken  from  Gotcha  data  set 
documentation  [25]. 


affect  the  suitability  of  the  data  for  use  in  the  following  empirical  studies,  as  long  as 
the  distance  from  the  antenna  phase  center  to  the  scene  center  is  known  for  every 
sample  in  azimuth. 

One  notable  benefit  of  the  parking  lot  scene  featured  in  this  data  set  is  that  it 
includes  fifteen  calibration  targets,  as  depicted  in  Figure  38  and  listed  in  Table  16. 
There  is  one  large  top  hat,  shown  in  Figure  39,  as  well  as  seven  trihedrals  and  seven 
dihedrals  of  varying  sizes  and  orientations.  The  calibration  targets  were  placed  in  a 
held  near  the  parking  lot  as  shown  in  Figure  40.  These  calibration  targets  are  useful 
for  evaluating  the  classification  accuracy  of  the  SPLIT  algorithm. 

Figure  41  presents  scatterer  classification  results  for  the  fifteen  calibration  targets 
of  the  Gotcha  data  set.  The  only  scatterers  that  appear  to  be  well  classified  are 
27TR1  near  (—7,52)  m,  DR3  near  (—18,33)  m,  DR5  near  (—13,32)  m,  and  DR7 
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Table  16.  Calibration  Targets  for  Gotcha  Data  Set  [25]. 


Calibration  Target 

ID# 

x  (m) 

y  (m) 

z  (m) 

15  in  Trihderal 

15TR-01 

-32.14 

42.54 

-0.53 

15  in  Trihderal 

15TR-03 

-28.09 

38.67 

-0.42 

15  in  Trihderal 

15TR-04 

-13.86 

37.70 

-0.05 

15  in  Trihderal 

15TR-05 

-24.39 

32.96 

-0.33 

15  in  Trihderal 

15TR-06 

-32.50 

33.41 

-0.57 

15  in  Trihderal 

15TR-07 

-5.12 

22.98 

-0.05 

27  in  Trihderal 

27TR-01 

-7.51 

51.47 

-0.09 

12  in  12  in  Dihedral 

DR-01 

-15.55 

42.96 

-0.13 

12  in  12  in  Dihedral 

DR-02 

-26.16 

45.64 

-0.43 

12  in  12  in  Dihedral 

DR-03 

-18.58 

33.53 

-0.18 

12  in  12  in  Dihedral 

DR-04 

-20.88 

27.10 

-0.23 

12  in  8  in  Dihedral 

DR-05 

-13.24 

32.09 

-0.09 

12  in  12  in  Dihedral 

DR-06 

-29.27 

24.48 

-0.48 

12  in  8  in  Dihedral 

DR-07 

-26.15 

17.50 

-0.44 

Figure  39.  Photo  of  top  hat  calibration  target  from  Gotcha  data  set.  [25]. 


near  (—26, 17)  m,  where  trihedral  and  dihedral  classifications  dominate.  Overall,  the 
classifier  showed  little  confusion  in  determining  the  polarimetric  response  of  scatterers, 
but  it  showed  some  confusion  in  determining  the  frequency  response  of  scatterers.  For 
instance,  none  of  the  trihedrals  is  classified  as  even-bounce,  and  none  of  the  dihedrals 
or  top  hat  is  classified  as  odd-bounce. 

For  this  data  set,  the  confusion  in  determining  frequency  response  is  most  likely 
due  to  a  combination  of  poor  coherency,  stray  clutter  energy,  unknown  system  biases, 
and  small  fractional  bandwidth.  Errors  in  measuring  the  distance  between  the  scene 
center  and  the  antenna  phase  center  disrupts  coherency  of  the  SAR  signal  and  vio- 
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Figure  40.  Photo  of  Gotcha  scene  [25].  The  calibration  targets  were  placed  inside  the 
area  marked  with  a  black  square. 


lates  the  assumptions  of  the  imaging  operator.  This,  in  turn,  violates  assumptions  of 
the  multi-peak  model  and  SPLIT  algorithm.  While  pass  number  six  showed  superior 
coherency  when  compared  to  the  other  passes,  there  still  appeared  to  be  significant 
errors  affecting  coherency,  even  after  using  the  autofocus  corrections  made  available 
in  the  data  set.  Furthermore,  the  parking  lot  scene  was  extracted  using  digital  spot¬ 
lighting  [135]  from  a  5  km  spotlight  SAR  collection,  and  it  appeared  that  significant 
energy  from  moving  vehicles  and  other  clutter  far  removed  from  the  parking  lot  had 
leaked  into  the  data.  Also,  it  is  unclear  how  best  to  normalize  the  data  for  possible 
system  induced  biases  in  the  data,  such  as  power  differences  due  to  the  antenna  gain 
pattern.  Finally,  fractional  bandwidth  of  6.5-percent  will  likely  lead  to  good  classi¬ 
fication  results  for  simple  calibration  targets  under  ideal  conditions.  However,  this 
fractional  bandwidth  was  far  too  narrow  to  overcome  the  significant  deficiencies  of 
poor  coherency,  stray  clutter  energy,  and  unknown  biases. 

These  non-ideal  conditions  provided  an  opportunity  to  stress  the  scatterer  classi¬ 
fier  in  way  that  the  previous  two  data  sets  did  not.  The  classification  results  for  the 
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Figure  41.  Scatterer  classification  results  for  the  Gotcha  calibration  targets.  Targets 
are  drawn  onto  the  image  using  white  lines,  where  the  trihedrals  and  dihedrals  are  at 
three  times  normal  scale  for  visibility.  The  top  40  dB  of  pixel  intensities  are  shown  in 
the  image. 


Gotcha,  data  set  demonstrate  some  of  the  limitations  of  the  SPLIT  algorithm  as  well 
as  the  concept  of  scatterer  classification  by  phase  history  decomposition,  in  general. 
Note  that  these  non-ideal  conditions  greatly  affect  classification  accuracy  using  image 
segmentation  methods,  and  no  satisfying  scatterer  classification  results  using  image 
segmentation  methods  have  been  published  for  the  Gotcha  Public  Release  Data  Set. 
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In  order  to  understand  how  bandwidth,  clutter,  and  interference  affect  classification 
accuracy,  the  following  section  presents  some  additional  analyses  and  experiments. 

5.5  Sensitivity  to  Bandwidth,  Clutter,  and  Interference 

To  further  illustrate  some  of  the  limitations  of  scatterer  classification  by  phase 
history  decomposition,  the  section  concludes  with  experiments  examining  the  effects 
of  bandwidth,  clutter,  and  interference  on  frequency  parameter  estimation.  Clas¬ 
sification  accuracies  for  ideal  scatterers  are  derived  from  Monte  Carlo  simulations 
based  on  a  signal  model  that  accounts  for  statistical  interference  from  neighboring 
canonical  scatterers,  clutter,  and  noise.  The  fractional  bandwidth  is  shown  to  be  a 
primary  contributor  to  classification  accuracy,  and  for  typical  SAR  applications,  frac¬ 
tional  bandwidths  of  10-percent  or  more  are  recommended  as  a  rule-of-thumb.  The 
following  studies  focus  on  the  frequency  parameter  a'  because  it  is  more  sensitive  to 
bandwidth  than  the  polarimetric  parameters.  The  assumption  is  that  signal  condi¬ 
tions  which  favor  a  good  a'  estimate  also  produce  excellent  polarimetric  parameter 
estimates. 


5.5.1  Signal  Model. 

SAR  imaging  is  a  coherent  process  that  can  be  interpreted  as  matched  filtering 
the  received  signal  to  the  expected  signal  due  an  ideal  point  scatterer  centered  at 
each  pixel  location  [76,  93].  Therefore,  by  replacing  the  integrals  in  Equation  (37)  of 
Section  4.4  with  summations  and  normalizing  the  integrand  by  HB^_f  )H0(e-e  ) ;  the 
subimage  peak  intensity  clue  to  a  single  canonical  point  scatterer  is  approximated  by 


fWn.  1  =  |S„(x„9,)|2  «  (iV,|^U/c)“'/2||/cl)2, 


(92) 
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where  the  narrow-band  approximation,  /  =  fc,  is  used,  and  Ns  is  the  number  of 
frequency  and  azimuth  samples.  For  wide-angle  and  wide-band  systems,  the  number 
of  samples  is  expected  to  be  quite  large.  For  instance,  in  the  case  of  128  frequency 
samples  and  128  azimuth  samples  the  total  number  of  samples  is  Ns  =  16384.  Note 
that  Equation  (92)  assumes  a  sampling  loss  of  zero,  whereby  the  location  of  the 
scatterer  is  exactly  in  the  center  of  a  pixel  [133].  The  zero  sampling  loss  assumption 
corresponds  with  use  of  the  SPLIT  algorithm,  whereby  canonical  scatterers  located  at 
or  near  the  center  of  a  pixel  have  a  greater  chance  of  passing  the  stable  peak  criterion. 

In  practice,  the  subimage  peak  intensity  also  includes  contributions  due  to  inter¬ 
ference  from  neighboring  scatterers,  clutter,  and  noise.  Therefore,  the  subimage  peak 
signal  model  is  given  as 


y  =  Ns 


A  (jfc)  //'  +  b interference  T  Sclutter  T  noise  \fc 


(93) 


where  the  interference,  clutter,  and  noise  are  modeled  as  random  variables  (RVs). 
These  are  discussed  next,  starting  with  the  noise  term. 


5.5. 1.1  Noise  Model. 

It  is  common  to  model  SAR  system  noise  as  additive  white  Gaussian  noise  (AWGN) 
[133,  93].  The  noise  has  a  circular  symmetric  complex  normal  distribution  expressed 
as  [133,  116] 

Snoise  =  w,  W  ~  CM{  0,  o\ y) ,  hd,  (94) 

where  samples  of  the  AWGN  are  independent  and  identically  distributed  (iid).  The 
magnitude  is  Rayleigh  distributed  with  variance,  VAR, (IT)  =  j^,  so  that  the 
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subimage  peak  intensity  clue  to  noise  alone  is 


incise  =  NsE  [|W|2]  |/c|2  =  Ns 


4  —  7T 


-a. 


w 


f2 

J  C  5 


(95) 


where  E[- ]  is  the  expected  value  of  the  RV. 

As  a  result,  the  signal-to-noise  ratio  (SNR)  is 

SNR  =  Psignal  «  N  2\A'\2f? 

P  ■  T)<72  ' 

1  noise  /l  )UW 


(96) 


It  can  be  seen  that  coherent  integration  serves  to  greatly  increase  the  SNR.  For 
instance,  in  the  case  of  large  Ns,  even  scattering  centers  with  very  weak  amplitudes 
can  obtain  a  high  SNR. 


5. 5. 1.2  Clutter  Model. 

Clutter  is  caused  by  backscatter  from  the  natural  environment  [133].  Examples 
include  surface  clutter  from  the  ground  and  volume  clutter  from  rain  [133].  The 
definition  of  clutter  depends  upon  the  unwanted  signal,  as  determined  by  the  specific 
application.  For  the  purposes  of  these  experiments,  the  clutter  is  modeled  as  surface 
clutter  from  a  large  region  of  unresolved  scatterers,  where  none  of  the  individual 
scatterers  is  significantly  stronger  than  the  others  [133]. 

Under  these  assumptions,  the  clutter  has  a  circular  symmetric  complex  normal 
distribution  expressed  as  [133,  116] 


flutter  =  \x I ej4,x,  \x\~  Rayleigh (o-x),  4>x~U[- 7r,7r],  X  1  <f>x,  (97) 


where  the  magnitude  has  a  Rayleigh  distribution  and  the  phase  difference  between  the 
signal  and  clutter,  <f>x,  is  iid  and  has  a  uniform  distribution.  Here,  the  magnitude  and 
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phase  are  independent  RVs.  Samples  of  clutter  are  identically  distributed,  but  unlike 
noise,  they  are  correlated.  Therefore,  the  subimage  peak  intensity  clue  to  clutter  alone 
is  greater  than  that  for  uncorrelated  samples 


-^clutter  NSE 


abs(|X|e^)2  |/c|2  =  N, 


7T 


12V2 


o 


x 


f, 


(98) 


and  less  than  that  for  perfect  correlation 


P flutter  <  N:E 


abs(|X|e^)2  | fc\2  =  Nt 


n 


I2V2 


o 


x 


/, 


(99) 


Here  the  factor  4^  produces  an  RMS  value  due  to  the  random  phase. 

Because  of  this  wide  range  of  possible  values,  the  clutter  power  is  parameterized, 
with  the  parameter  being  measured  or  approximated  for  a  given  type  of  surface.  Fol¬ 
lowing  Reference  [133],  the  parameterized  clutter  model  assumes  perfect  correlation, 
and  then  simply  scales  the  power  using  an  effective  clutter  RCS,  oc,  for  the  clutter 
occupying  the  area,  Ac,  illuminated  by  a  single  range  cell  of  the  radar.  These  are  com¬ 
bined  to  provide  a  normalized  clutter  RCS,  o°  =  with  values  generally  ranging 
from  -10  dB  for  mountains  and  urban  terrain  to  -40  dB  for  grassland,  depending  upon 
depression  angle  and  operating  frequency  [133].  The  signal-to-clutter  ratio  (SCR)  is 
[133], 


SCR  = 


Signal  „  AC|Hf  /c“72 


Or 


P 

1  clutter 

where  using  Equation  (99),  the  relation  to  the  clutter  RV  is 


(100) 


=  Aco°  =  N 
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(101) 
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5.5. 1.3  Interference  Model. 


For  these  experiments,  interference  is  caused  by  the  sidelobes  of  the  PSFs  of 
strong  neighboring  scatterers  in  the  scene.  These  sidelobes  are  coherently  summed  to 
produce  the  interference  term  for  a  given  sample.  The  difference  between  interference 
and  clutter  is  that  the  clutter  model  assumes  none  of  the  individual  scatterers  are 
significantly  stronger  than  the  others,  while  the  interference  model  does  not  assume 
uniformity  among  scatterers  [133,  116]. 

Under  these  assumptions,  the  interference  is  modeled  as  having  a  magnitude  with 
a  gamma  distribution  and  a  phase,  (f>z,  with  a  uniform  distribution.  The  interference 
is  modeled  as 


interference  =  \Z\ e34>z\  \  Z\  ~  T(Ni,  <7Z) ,  <Pz  ~  W[-7T,  7r] ,  Z  X  <j>Z.  (102) 


According  to  Reference  [116],  the  gamma  distribution  parameters  can  be  interpreted 
as  Ni  being  equal  to  the  number  of  dominant  interfering  scatterers  with  mean  RCS 
equal  to  az.  The  samples  of  interference  are  identically  distributed  and,  for  simplicity, 
assumed  to  be  perfectly  correlated  in  magnitude  and  iid  in  phase.  Therefore,  the 
subimage  peak  intensity  due  to  interference  alone  is 


-^interference  ^  Ah  E 


abs(|Z|e^)2  |/c|2  =  N\ 
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L  V2 


f2 

J  C  1 


where  the  factor  produces  an  RMS  value  due  to  the  random  phase. 
Thus,  the  corresponding  signal-to-interference  ratio  (SIR)  is 


SIR 


Pi 


interference 


Ni(r2z 


(103) 


(104) 
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5.5.2  Scatterer  Classification  Accuracy  Experiments. 

The  curve  fitting  algorithm  described  in  Section  5. 1.3.1  allows  the  frequency  pa¬ 
rameter  to  assume  a  value  in  the  continuous  interval  a'  G  [—4,4],  This  was  done 
to  facilitate  weighted  averaging  of  a'  from  multiple  subapertures  and  co-polarization 
channels.  In  contrast,  the  experiments  presented  in  this  section  are  restricted  to  a 
single  aperture,  where  the  correlation  is  unspecified  between  subapertures  and  co¬ 
polarization  channels  for  the  signal  model  in  Equation  (93).  In  this  case,  it  is  simpler 
to  allow  the  frequency  parameter  to  assume  a  discrete  value,  a'p  =  p ,  in  the  set 
p  G  {—5,  —3,  —2,  —1,  0, 1,  2,  3, 5}.  In  this  way,  the  classification  decision  is  made  ac¬ 
cording  to  the  minimum  total  least  squares  metric  between  the  simulated  observation 
vector  and  the  seven  ideal  curves  given  by  f  (p)  from  Equation  (62)  of  Section  5. 1.3.1, 
for  a  given  fractional  bandwidth,  Q  =  In  this  case,  a  classification  decision  at 
the  extremes  of  p  G  {—5,5}  causes  the  SPLIT  algorithm  to  erroneously  reject  the 
canonical  scatterer  as  non-canonical,  that  is  a'  (f  [—4,4],  Such  an  event  can  be  in¬ 
terpreted  as  a  missed  detection  due  to  a  high  level  of  interference,  clutter,  noise,  or  a 
combination  of  these  in  the  simulated  observation  vector. 

Referring  to  Equation  (63)  of  Section  5. 1.3.1  and  to  Equation  (93),  the  mth  sim¬ 
ulated  observation  vector  for  a  given  a',  SIR,  SCR,  and  SNR  is  given  by 

<Tm  =  [\y\lm,\y\lm,---,\y\lm]T 

(|A'|/iW/2  +  Zme?**  i/cl  +  Xme?+xifcl  +  W[™lfc  i) 

(\A'\f£a'/2  +  zme^fc2  +  X.me^fc2  +  W™lfc2) 

(|^|/c/+a72  +  Zme^fcI  +  Xme^fcI  +  W^hi) 

where  Zm  and  Xrn  are  the  values  of  the  RVs  Z  and  X,  respectively,  for  the  mth 
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observation.  Likewise,  (f>zt  and  4>xt  are  the  values  of  the  RVs  <pz  and  (f>x  for  the  ith 
subband.  Finally,  the  real  component  of  the  noise  is  taken  to  be  in-phase  with  the 
signal,  without  loss  of  generality,  so  that  W[^al  is  the  value  of  the  in-phase  component 
of  the  RV  W  for  the  ith  subband  and  mth  observation.  Note  that  the  normalization 
factor,  um,  includes  the  factors  ./V2/2,  which  assumes  that  the  number  of  phase  history 
samples  is  identically  equal  to  Ns  for  each  subdomain. 

The  experiments  in  this  section  normalize  \A'\  and  fc  to  unity,  without  loss  of 
generality.  Also,  to  normalize  the  signal  power,  the  frequency  parameter  is  set  to 
a!  =  0,  which  represents  an  average  value  assuming  equal  prior  probabilities  for  each 
class  of  scatterer.  As  a  result,  the  variances  are  normalized  to 

__  2  s/2  _  2  2\[2  _  2  2  Ns  (  ^ 

^  =  iV-(SlR)  ’  °x=  (4-tt)(SCR)’  ^  (4-tt)(SNR)-  (106) 

Because  AWGN  is  uncorrelated,  the  SNR  is  usually  much  larger  than  SIR  or  SCR. 
Therefore,  it  is  ignored  in  these  experiments  by  setting  SNR  =  oo.  This  allows  the 
primary  effects  of  bandwidth,  SCR,  and  SIR  to  be  examined  in  more  detail.  The 
classification  accuracy  for  varying  fractional  bandwidths,  SIR,  and  SCRo  are  shown 
in  Figure  42.  These  were  produced  from  Monte  Carlo  simulations  of  ten  thousand 
observation  vectors  for  each  combination  of  fractional  bandwidth,  SIR,  and  SCRo- 
The  use  of  SCRo  accounts  for  the  fact  that  SCR  is  actually  dependent  upon  subimage 
pixel  area,  Ac,  which  in  turn  is  dependent  upon  bandwidth.  In  order  normalize  the 
SCR  for  each  fractional  bandwidth,  SCRo  was  chosen  according  to  the  Gotcha  data 
set,  where  Ac  =  0.45  m2  and  flo  =  0.065.  Hence,  in  order  to  account  for  changing 
subimage  pixel  area,  the  SCR  for  each  fractional  bandwidth  varied  as 

SCR(j3)  =  SCR0(/5o//32).  (107) 
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Figure  42.  Classification  accuracies  for  the  case  of  I  =  3  subbands,  (3q  =  0.065,  Ni  =  2 
interferers,  and  SNR  =  oo  for  varying  fractional  bandwidth,  (3,  SIR,  and  SCRq. 
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In  addition,  the  number  of  subbands  was  set  to  /  =  3,  and  the  number  of  neighboring 
scatterers  was  set  to  iVj  =  2  for  these  experiments. 

Correct  classifications  occurred  when  the  declared  a'  was  equal  to  the  true  a' , 
that  is  when  a!  =  p.  Alternately,  note  that  the  canonical  range,  a'  G  [—4,4],  means 
that  a  declared  value  of  cd  =  —3  for  true  a'  —  —  2  and  a  declared  value  of  a'  =  3 
for  true  a'  =  2  could  both  be  recategorized  as  correct  classifications,  if  desired.  In 
these  cases,  the  extracted  feature  vector  contributes  to  the  average  feature  vector  so 
as  to  produce  a  correct  classification,  but  with  slightly  lower  measure  of  fitness.  In 
this  alternate  case,  the  classification  accuracies  in  Figure  42  would  increase. 

Figure  42  reveals  that  classification  accuracy  is  greatly  impacted  by  fractional 
bandwidth,  regardless  of  SIR  or  SCR.  For  instance,  Figure  42(a)  reveals  that  for 
very  low  fractional  bandwidth,  (3  =  0.02,  SIR  and  SCR0  must  both  be  at  60  dB 
or  greater  before  the  classification  results  are  even  above  25-percent.  Any  point 
where  classification  accuracy  is  above  20-percent  can  be  loosely  interpreted  as  a  case 
where  the  least  squares  classifier  produces  better  results  than  a  random  classifier 
resembling  to  a  ‘coin  toss.’  However,  the  inclusion  of  p  G  {—5,  —3,3,5}  complicates 
the  interpretation. 

Note,  that  the  Euclidean  norm  varies  slightly  between  the  adjacent  pairs  of  ideal 
curves,  ||f(p)  —  f (p  +  1)||2.  For  example,  in  the  case  of  /  =  3  and  j3  =  0.5, 

||f(— 2)  -  f(-l)||2  =  0.1768, 

||f(— 1)  —  f(0)||2  =  0.1782, 

(108) 

||f(0)  -  f(l)||2  =  0.1849, 

1 1  f  (1)  —  f  (2)  1 1 2  =  0.1967. 

This  is  one  factor  that  will  cause  classification  accuracy  to  vary  slightly  for  different 
values  of  a' .  These  can  be  adjusted  to  become  more  equal  by  adjusting  the  normal- 
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ization  factor,  fc,  used  in  (/)  =  p.  With  a  small  adjustment,  it  is  possible  that  overall 
classification  accuracy  can  be  slightly  improved  for  these  experiments.  However,  the 
signal  model  includes  assumptions,  such  as  non-dispersive  clutter  and  interference, 
which  reduce  the  value  of  such  fine  precision  adjustments,  in  practice.  Furthermore, 
the  optimal  normalization  frequency  is  expected  to  be  closer  to  fc,  in  general,  for 
rectangularly  shaped  phase  histories  where  the  Jacobian,  \fc\,  is  removed  from  the 
imaging  operator,  as  discussed  in  Section  5.1.3. 

The  confusion  matrices  in  Table  17  are  provided  to  illustrate  how  classification 
accuracy  varies  between  different  values  of  a'  for  the  cases  of  interference  only  (top) 
and  clutter  only  (bottom).  The  values  in  red  italics  indicate  the  number  of  correctly 
classified  scatterers  out  of  one-hundred  thousand  trials  for  each  value  of  a'.  There 
seems  to  be  competing  trends  in  the  data.  While  in  some  cases  classification  accuracy 
tends  to  improve  as  a'  increases,  in  other  cases,  accuracy  improves  as  a'  decreases. 
The  inflexion  points  for  these  trends  depend  upon  the  amount  of  clutter  or  interference 
present.  In  theory,  these  trends  can  be  modified  by  adjusting  the  normalization  factor 
in  (/)  =  p,  but  since  the  SCR  and  SIR  are  not  usually  known  a  priori ,  this  is  not 
very  practical. 

The  confusion  matrix  in  the  bottom  left-hand  corner  of  Table  17  for  SCR  =  30  dB 
seems  to  mimic  the  marginal  classification  accuracy  for  the  Gotcha  data  experiment 
of  Figure  41  of  Section  5.4.4.  However,  the  SCR  for  the  top  hat  and  six  smaller 
trihedrals  can  be  readily  calculated  as 


SCR 


(Ttophat  _  800  m2 

a°Ac  ~  (-10  dB)  (0.455  m2) 


42  dB, 


(109) 


where  atophat  ~  800  m2  is  approximated  using  the  equation  for  cylinder  RCS  given 
in  Figure  8  of  Section  2.2.1  and  a0  ~  —10  dB  is  taken  from  table  7.11  of  Reference 
[133].  The  difference  between  an  SCR  of  42  dB  and  a  classification  performance  that 
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Table  17.  Simulated  Confusion  Matrices  for  the  Case  of  (3o  =  0.065,  I  =  3  subbands, 
Ni  =  2  interferers,  and  SNR  =  oo. 


True  a!  value 

SIR  =  30  dB,  SCR  =  oo 

SIR  =  50  dB,  SCR  =  oo 

-2 

-1 

0 

1 

2 

-2 

-1 

0 

1 

2 

Declared  a!  value 

-5 

21883 

15463 

11130 

7906 

5863 

64 

4 

- 

1 

- 

-3 

16654 

10607 

7210 

5202 

3706 

8067 

334 

18 

- 

- 

-2 

22886 

12389 

7808 

5374 

3612 

83929 

7904 

328 

20 

1 

-1 

12470 

22764 

12540 

7832 

5325 

7579 

83722 

7728 

372 

26 

0 

7632 

12467 

22771 

12500 

7845 

343 

7701 

83729 

7792 

327 

1 

5301 

7804 

12333 

22517 

12593 

18 

322 

7857 

83604 

7723 

2 

3707 

5341 

7775 

12545 

22682 

- 

12 

323 

7859 

83713 

3 

3606 

5127 

7259 

10754 

16847 

- 

1 

17 

350 

8132 

5 

5861 

8038 

11174 

15370 

21527 

- 

- 

- 

2 

78 

Classification  Accuracy  =  22.7240% 

Classification  Accuracy  =  83.7394% 

True  a'  value 

SIR  =  oo,  SCR  =  30  dB 

SIR  =  oo,  SCR  =  50  dB 

-2 

-1 

0 

1 

2 

-2 

-1 

0 

1 

2 

Declared  a!  value 

-5 

28514 

22183 

17310 

13022 

10194 

59 

1 

- 

- 

- 

-3 

13599 

10041 

7797 

6455 

5059 

13213 

499 

2 

- 

- 

-2 

15940 

10052 

7190 

5630 

4454 

73425 

12723 

529 

6 

- 

-1 

9965 

15796 

9997 

7162 

5624 

12789 

73500 

12915 

525 

6 

0 

7151 

9875 

15787 

10195 

7073 

512 

12769 

73144 

12824 

496 

1 

5492 

7049 

9891 

15766 

10320 

2 

502 

12899 

73233 

12774 

2 

4323 

5570 

7237 

9991 

15566 

- 

6 

508 

12863 

73262 

3 

5049 

6344 

7785 

10072 

13762 

- 

- 

3 

549 

13404 

5 

9967 

13090 

17006 

21707 

27948 

- 

- 

- 

- 

58 

Classification  Accuracy  =  15.771% 

Classification  Accuracy  =  73.3128% 

mimics  an  SCR  of  30  dB  indicates  that  there  is  approximately  a  12  dB  loss  factor 
in  the  Gotcha  data  due  to  poor  coherency  and  stray  clutter,  as  discussed  earlier  in 
Section  5.4.4. 


5.5.3  Coupling  Between  Bandwidth  and  Interference. 

As  a  simplification,  the  signal  model  and  experiments  treat  bandwidth  and  inter¬ 
ference  as  independent  factors  affecting  classification  accuracy.  However,  in  reality, 
these  are  actually  interdependent  factors.  The  amount  of  sidelobe  interference  af¬ 
fecting  the  subimage  pixel  peak  intensity  is  dependent  upon  the  size  of  the  pixel. 
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This  coupling  of  interference  power  and  bandwidth  was  not  explicitly  accounted  for 
earlier,  but  can  have  a  significant  impact  on  classification  accuracy,  as  illustrated  in 
the  following  experiments. 

Figure  43  illustrates  the  impact  of  bandwidth-interference  coupling  on  classifica¬ 
tion  accuracy  for  D2.  All  data  parameters  from  Table  13,  including  center  frequency, 
are  held  constant  except  for  fractional  bandwidth.  The  main  source  of  interference 
is  due  to  target-pedestal  interactions  in  the  center  of  the  image.  As  the  fractional 
bandwidth  decreases,  the  SIR  decreases  for  pixels  at  and  near  the  center  of  the  image. 
As  a  result,  the  classification  accuracy  for  these  pixels  suffers.  A  most  obvious  exam¬ 
ple  appears  in  Figure  43(d),  where  the  ring  of  tri/dihedral  classifications  from  Figure 
43(e)  disappears.  Because  this  object  is  so  simple  and  the  edge  scatterers  are  rela¬ 
tively  isolated  from  neighboring  scatterers,  the  edge  classifications  remain  accurate 
even  for  smaller  fractional  bandwidths. 

Figure  44  illustrates  the  impact  of  bandwidth-interference  coupling  for  the  Toyota 
Tacoma.  The  Toyota  Tacoma  is  a  more  complex  target  than  D2,  causing  a  lower 
SIR  for  many  of  the  scatterers.  As  a  result,  the  major  features  described  in  Section 
5.4.3  are  evident  in  Figure  44(c),  but  dissipate  at  smaller  fractional  bandwidths. 

Note  that  the  SCR  and  SNR  are  approximately  infinite  in  Figures  43  and  44, 
so  that  the  impact  of  interference  could  be  examined  in  isolation.  Based  on  these 
and  previous  examples,  it  is  recommended  that  the  SPLIT  algorithm  be  used  only 
for  fractional  bandwidths  of  10-percent  or  more  for  typical  SAR  applications,  where 
targets  are  comprised  of  multiple  canonical  scatterers  in  close  proximity. 

5.5.4  Multiple  Observations  and  Oversampling. 

Classification  accuracy  is  expected  to  improve  as  multiple  observations  of  the  scat- 
terer  are  available.  Recall  that  multiple  observations  are  combined  through  a  weighted 
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Figure  43.  Classification  results  for  D2  with  changing  fractional  bandwidth. 
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Figure  44.  Classification  results  for  the  Toyota  Tacoma  with  changing  fractional  band¬ 
width. 


average  of  feature  vectors,  as  described  in  Section  5.2.  For  the  frequency  parameter, 
multiple  observations  are  potentially  available  through  the  two  co-polarization  chan¬ 
nels,  HH  and  VV,  and  through  multiple  subapertures.  For  co-polarization  channels, 
the  observations  are  considered  completely  independent,  while  for  multiple  subaper¬ 
tures,  the  observations  are  correlated  because  the  clutter  and  interference  contribu¬ 
tions  are  correlated.  For  two  independent  observation  vectors  of  equal  weight,  the 
benefit  to  classification  accuracy  can  be  equated  to  a  doubling  of  the  SIR  and  SCR 
for  the  single  observation  case.  Although  multiple  observations  are  common  in  prac¬ 
tice,  these  usually  do  not  produce  independent  observation  vectors  of  equal  weight. 
Therefore,  improvements  to  classification  accuracy  can  be  expected  when  multiple 
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observations  are  available,  but  usually  less  than  that  for  the  ideal  scenario  of  dou¬ 
bling  of  the  effective  SIR  and  effective  SCR  as  the  number  of  observations  doubles. 
The  confusion  matrices  in  Section  5.5.2  considered  only  a  single  observation,  while 
the  annotated  images  in  Section  5.5.3  include  all  available  opportunities  for  multiple 
observations. 

Another  way  to  increase  the  number  of  observations  is  by  decomposing  the  phase 
history  into  an  increasing  number  of  subbands.  In  this  case,  the  clutter  and  inter¬ 
ference  is  expected  to  be  highly  correlated  among  the  subimages,  but  the  effective 
SIR  and  effective  SCR  are  expected  to  improve  slightly  none  the  less.  For  example, 
Figure  45  shows  the  classification  accuracy  for  the  case  of  /  =  5  and  reveals  that 
classification  accuracy  only  slightly  improves  compared  to  the  case  of  /  =  3  in  Figure 
42.  The  model  used  in  this  experiment  does  not  include  the  expected  increase  in  cor¬ 
relation  between  interference  and  clutter  samples  when  I  =  5;  so  any  improvements 
in  classification  accuracy  represent  a  best  case  scenario. 

The  last  consideration  for  classification  accuracy  is  the  effect  of  pixel  oversampling. 
Pixel  oversampling  is  defined  as  increasing  the  number  of  pixels  in  each  resolution  cell, 
and  can  be  interpreted  as  an  interpolation  of  the  image  samples.  Because  the  SPLIT 
algorithm  extracts  feature  vectors  on  a  pixcl-by-pixcl  basis,  the  size  and  number  of 
the  pixels  could  potentially  have  an  impact  on  classification  accuracy.  The  following 
is  only  a  cursory  discussion  on  classification  accuracy  trends  due  to  oversampling 
because  it  is  difficult  to  express  the  impact  of  oversampling  analytically.  Furthermore, 
the  use  of  oversampling  is  often  limited  in  practice  because  it  results  in  an  exponential 
increase  in  processing  costs. 

Oversampling  limits  the  amount  of  interference  energy  in  a  given  pixel  and  is 
therefore  expected  to  increase  the  effective  SIR.  However,  with  oversampling  there 
are  more  pixel  boundaries  and  greater  opportunity  for  a  scatterer  to  be  located  very 
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Figure  45.  Classification  accuracies  for  the  case  of  /  =  5  subbands,  (3q  =  0.065,  Ni  =  2 
interferers,  and  SNR  =  oo  for  varying  fractional  bandwidth,  (3,  SIR,  and  SCRq. 
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near  a  pixel  boundary.  When  a  scatterer  is  located  very  near  a  pixel  boundary, 
there  is  a  chance  that  noise  and  clutter  interference  may  cause  the  subimage  peak 
to  migrate  into  a  neighboring  pixel,  thereby  causing  the  scatterer  to  fail  the  stable 
peak  criterion  of  the  SPLIT  algorithm,  as  discussed  in  Section  5.1.1.  Therefore, 
oversampling  increases  the  probability  of  missed  detections,  which  in  turn,  decreases 
the  number  of  observations.  Alternately,  for  scatterers  near  a  boundary,  the  sampling 
loss  may  vary  for  each  subband,  thereby  skewing  the  estimation  of  the  frequency  and 
polarimetric  parameters  during  feature  vector  extraction.  As  a  result,  experience 
shows  that  classification  accuracy  may  show  only  a  slight  overall  improvement  with 
increased  oversampling,  but  usually  not  enough  to  warrant  the  increased  processing 
cost. 

As  an  example,  classification  results  for  2x  oversampling  are  given  in  Figure  46 
for  D2.  When  these  classification  results  are  compared  to  those  in  Figure  43,  there 
does  not  seem  to  be  a  significant  increase  in  classification  accuracy  due  to  oversam¬ 
pling.  Furthermore,  missed  detections  due  to  oversampling  are  evident,  particularly 
in  comparing  Figure  46(e)  to  Figure  43(e). 
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Figure  46.  Classification  results  for  D2  with  2x  oversampling  and  changing  fractional 
bandwidth. 
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VI.  Integrated  Domain  Decomposition  and  Scatterer 
Classification  Algorithm 


This  chapter  explains  how  the  subimages  used  for  scatterer  classification  can  be 
efficiently  combined  to  approximate  a  conventional  SAR  image  reconstructed  from  the 
original  phase  history.  Alternately,  it  is  always  possible  to  perform  SPLIT-based  scat¬ 
terer  classification  separate  from  conventional  SAR  imaging,  but  this  is  less  efficient. 
Therefore,  it  is  of  interest  to  reconstruct  the  SAR  image  directly  from  the  subimages, 
using  domain  decomposition  imaging  techniques.  This  integrated  approach  increases 
the  overall  computational  efficiency  and  usefulness  of  SPLIT-based  scatterer  classifi¬ 
cation.  Although  domain  decomposition  imaging  introduces  some  imaging  error  due 
to  interpolation  and  subwindow  summation,  the  error  is  controllable. 

The  discussion  and  evaluation  of  the  integrated  algorithm  is  organized  as  follows. 
Section  6.1  explains  how  the  coarse  resolution  subimages  are  be  interpolated  and 
combined  to  approximate  a  conventional,  fine  resolution  image.  The  imaging  error 
due  to  interpolation  is  shown  to  be  controllable,  and  example  SAR  images  illustrate 
how  the  order  of  the  imaging  operator  is  expected  to  affect  imaging  accuracy.  Section 
6.2  presents  subwindow  design  principles  for  approximating  full  band  and  full  aperture 
windows  used  in  conventional  imaging.  Additional  analysis  reveals  that  the  imaging 
error  clue  to  the  approximation  is  especially  well-controlled  for  full  360°  apertures, 
and  I  >  5  subbands.  However,  for  other  cases,  some  imaging  artifacts  may  be  visible. 
Section  6.3  introduces  SAR  surveillance  applications  and  the  idea  of  coverage  area  as 
a  performance  metric.  It  presents  analysis  of  the  computational  cost  of  the  integrated 
algorithm  and  concludes  with  an  efficiency  study  using  the  Gotcha  data  set.  As  a 
result,  surveillance  SAR  applications  are  given  specific  consideration  throughout  this 
chapter. 
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Figure  47  provides  an  overview  of  the  integrated  algorithm.  Starting  in  the  upper 
left-hand  corner,  the  subaperture  windows  divide  the  phase  history  into  J  uniformly- 
spaced,  overlapping  subapertures.  Next,  in  the  upper  right-hand  corner,  the  halfband 
windows  divide  each  subaperture  into  I  =  3  uniformly-spaced,  overlapping  subbands. 
In  the  largest,  dotted-line  box  three  coarse-resolution  subimages  are  reconstructed 
from  the  subdomains  in  each  subaperture.  Two  processing  paths  proceed  downward 
from  the  subimages.  To  the  extreme  left,  the  SPLIT  algorithm  consists  of  feature 
vector  array  extraction,  summation,  and  classification,  as  described  in  Chapter  V.  In 
the  center,  domain  decomposition  imaging  consists  of  an  interpolation  and  weighted 
sum  of  the  subimages  to  form  a  subaperture  image.  Then,  the  subaperture  images 
are  summed  to  form  the  final  image.  Finally,  in  the  bottom  left-hand  corner,  the 
image  is  annotated  with  scatterer  classification  results.  The  steps  in  the  domain 
decomposition  imaging  path  are  discussed  in  more  detail  below. 

6.1  Subimage  Interpolation 

This  section  presents  a  detailed  interpolation  operator  for  the  integrated  algo¬ 
rithm.  The  interpolation  operator  assumes  overlapping  subimages  in  accordance  with 
a  SPLIT-based  decomposition  of  the  phase  history,  as  described  in  Chapter  V.  Refer¬ 
ring  to  Equation  3  of  Section  2. 1.2.2,  a  discrete  version  of  the  interpolation  operator 
is  given  as 

fey)  =  ^2cj^2ciI  {gij{x'  ,  (no) 

j  i 

where  Z{-}  is  the  interpolation  operator  and  the  ^  symbol  indicates  that  the  result 
is  an  approximation  to  the  fine  resolution  image  reconstructed  conventionally  from 
the  full  domain  of  the  phase  history.  Here,  the  subimage  pixel  coordinates  are  (x',  y') 
and  the  conventional  image  pixel  coordinates  are  (x, y).  Note  that  in  this  section, 
images  are  no  longer  expressed  as  continous  functions. 
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Figure  47.  A  detailed  diagram  of  the  integrated  scatterer  classification  and  imaging 
algorithm.  The  diagram  allows  visualization  of  the  effects  of  subdomain  windowing  for 
a  full  360°  aperture  and  I  =  3  subbands. 
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One  complication  of  SAR  image  interpolation  is  that  SAR  images  exhibit  rapid 
phase  variations  between  pixels  due  to  the  fact  that  SAR  imagery  is  recovered  from 
band  limited,  frequency  offset  data  [113].  Hence,  the  interpolation  operator  must 
remove  the  rapid  phase  variations  between  subimage  pixels  in  order  to  promote  accu¬ 
rate  interpolation  [12].  Recall  that  the  radar  system  collects  samples  of  the  scattered 
electric  held  at  different  aspect  angles  and  frequencies  as  it  moves  along  its  flight  path. 
For  any  single  point  scatterer  in  the  surveillance  area,  the  returned  phase  history  is 
exp (—j2kd(9)),  where  k  =  2irf  /c  is  the  wavenumber,  c  is  the  speed  of  light,  and  d{6) 
is  the  distance  to  the  airborne  radar  for  a  given  aspect  angle,  6  [76].  Considering 
that  each  image  pixel  corresponds  to  a  point  in  the  surveillance  area,  the  grid  of  pixel 
locations  produces  an  array  of  distances,  dxy.  Thus,  for  a  given  subdomain,  an  array 
of  central  phase  factors  is  given  by  A ^  =  exp(—  j2kcidx>y'(0cj)),  where  kcl  =  2iTfci/c 
is  the  central  wavenumber  of  the  ith  subband  and  d x'y'(0Cj)  is  an  array  of  distances 
between  the  airborne  radar  and  the  pixels  located  at  positions  (x',y'),  as  measured 
at  the  central  angle,  9Cj,  of  the  jth  subaperture. 

Following  [12],  the  interpolation  operator  multiplies  each  subimage  pixel  array  by 
its  corresponding  array  of  central  phase  factors  to  smooth-out  the  phase  variations 
between  pixels  before  2D  interpolation.  After  2D  interpolation  is  complete,  the  rapid 
phase  variations  are  reestablished  by  multiplying  by  the  conjugate  of  the  array  of 
central  phase  factors,  A*^-  =  exp(j2/cCjdxy(6)C:;)),  where  (x,  y)  are  the  positions  of  the 
interpolated  pixels.  Thus,  the  resultant  produced  by  the  interpolation  operator,  X{-}, 
is  succinctly  expressed  as 


y)  =  Z'{&i(x,,y');x,y}  =  a^ia^^x',/),  (in) 

where  I  performs  a  2D  interpolation  which  effectively  increases  the  sampling  rate 
in  both  coordinates  of  the  spatial  domain.  For  example,  recalling  that  the  SPLIT 
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Figure  48.  SAR  imagery  for  multilook  conventional  CBP  imaging  (left)  and  domain 
decomposition  CBP  imaging  with  a  subaperture  width  0  =  2°.  The  order  of  the  in¬ 
terpolation  increases  to  the  right.  It  appears  that  cubic  interpolation  is  required  to 
reduce  image  artifacts  to  an  indiscernable  level  compared  to  the  conventional  image. 


algorithm  uses  halfband  subwindows  in  frequency  and  assuming  that  the  full  aperture 
is  at  least  twice  the  size  of  the  subapertures;  the  interpolation  operator  effectively 
doubles  the  sampling  rate  in  each  dimension  so  that  the  resultant  has  four  times 
more  pixels  than  the  subimages.  The  ^  symbol  in  Equation  (111)  indicates  that  the 
interpolation  operator  introduces  a  controllable  error  [12], 


6.1.1  Imaging  Error  Due  to  Interpolation. 

The  imaging  error  due  to  interpolation  is  controllable  by  increasing  the  amount  of 
oversampling,  increasing  of  the  order  of  the  interpolator,  or  both  [12,  7].  Experience 
shows  that  for  very  narrow  subapertures,  linear  interpolation  is  both  efficient  and 
sufficient  to  control  interpolation  error  in  the  integrated  algorithm  [12],  However,  for 
larger  subapertures,  higher  order  interpolators  may  be  required.  For  example,  Figure 
48  shows  the  results  of  multilook  images  for  conventional  imaging  and  domain  de¬ 
composition  imaging  with  varying  orders  of  interpolation.  By  comparing  the  domain 
decomposition  images  to  the  conventional  image  on  the  left,  it  is  possible  to  visualize 
image  artifacts.  The  results  show  that  for  a  2°  subaperture,  cubic  interpolation  is 
required  to  reduce  image  artifacts  to  a  visually  indiscernible  level.  In  addition,  Figure 
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Figure  49.  SAR  imagery  for  multilook  conventional  CBP  imaging  (left)  and  domain 
decomposition  CBP  imaging  with  a  subaperture  width  0  =  10°.  It  appears  that  ideal 
interpolation  is  required  to  reduce  image  artifacts  to  an  indiscernable  level. 


(a)  Conventional  (b)  Nearest  (c)  Linear  Interp.  (d)  Cubic  Interp.  (e)  Ideal  Interp. 

Image.  Neighbor 

Interp. 


Figure  50.  SAR  imagery  for  conventional  CBP  imaging  (left)  and  domain  decompo¬ 
sition  CBP  imaging  with  a  subaperture  width  0  =  2°.  It  appears  that  even  nearest 
neighbor  interpolation  may  be  sufficient  for  some  applications. 


49  reveals  that  for  a  10°  subaperture,  ideal  interpolation  is  required.  In  this  case, 
ideal  interpolation  for  band  limited  signals  is  performed  by  taking  the  IFFT  of  the 
zero-padded  resultant  of  an  FFT. 

It  turns  out  that  coherent  summation  imaging,  as  expressed  in  Equation  (4)  of 
Section  2. 1.2. 4,  may  be  more  tolerant  to  interpolation  error  than  multilook  imaging, 
as  expressed  in  Equation  (5)  of  Section  2. 1.2.4.  For  example,  Figures  50  and  51  reveal 
that  even  nearest  neighbor  interpolation  may  be  sufficient  for  some  applications. 
Note  that  all  domain  decomposition  images  in  this  section  used  1  =  5  subbands. 
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(a)  Conventional  (b)  Nearest  (c)  Linear  Interp.  (d)  Cubic  Interp.  (e)  Ideal  Interp. 

Image.  Neighbor 

Interp. 


Figure  51.  SAR  imagery  for  conventional  CBP  imaging  (left)  and  domain  decompo¬ 
sition  CBP  imaging  with  a  subaperture  width  0  =  10°.  It  appears  that  even  nearest 
neighbor  interpolation  may  be  sufficient  for  some  applications. 


6.2  Subwindow  Design  and  Weighting 


In  addition,  tapered  subwindows  in  azimuth  improve  the  sidelobe  conditioning  of 
subaperture  images  used  in  video  SAR,  while  tapered  subwindows  in  both  azimuth 
and  frequency  provide  sidelobe  conditioning  to  improve  SIR  for  scatterer  classification. 
However,  for  the  integrated  algorithm,  these  subwindows  cause  artifacts  in  the  final 
image  when  their  summations  do  not  well- approximate  the  desired  full  domain  win¬ 
dows.  The  effects  of  data  windowing  can  be  visualized  in  the  detailed  diagram  of  the 
integrated  algorithm  provided  in  Figure  47.  Fortunately,  the  error  due  to  subdomain 
windowing  is  controllable  through  careful  design  and  use  of  Hanning  subwindows. 

The  following  discussion  addresses  the  controllable  error  in  each  spectral  dimen¬ 
sion,  separately.  The  approximation  is  expressed  as  a  summation  of  subwindows 

M 

*(£  -  tc)  ~  +  & m ),  (112) 

m= 1 

where  H<&  (£  —  £c)  is  either  the  full  aperture  window  in  azimuth  or  the  fullband  win¬ 
dow  in  frequency  with  region  of  support  <f>.  In  this  case,  M  6  N  is  the  number  of 
subwindows  with  region  of  support,  C  <3>,  in  that  dimension.  The  subwindows  are 
weighted,  scaled,  and  shifted  by  cm  G  M,  fi  G  M,  n  >  1,  and  £cm  G  M,  respectively. 
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The  window  error  for  Equation  (112)  is  defined  using  the  Euclidean  norm, 

M 

c-  ^  (us) 

m= 1 

and  is  a  function  of  the  type,  number,  spacing,  and  weighting  of  the  subwindows. 
For  simplicity  of  notation,  the  mth  shifted  window  is  H^m  =  +  £cm)  from 

Equation  (112).  Because  the  spectral  and  spatial  domains  form  a  Hilbert  space,  a 
solution  which  minimizes  the  window  error  in  the  spectral  domain  will  likewise  min¬ 
imize  artifacts  in  the  image  [15].  Raised  cosine  family  of  window  functions,  such  as 
Hanning,  Hamming,  and  Taylor,  are  often  used  in  SAR  imaging  [76],  but  other  sym¬ 
metric,  tapered  windows,  such  as  those  described  in  the  seminal  paper  by  Harris  [67], 
could  also  be  used.  Because  of  the  summation  in  Equation  (113),  the  window  error 
is  best  controlled  when  the  subwindows  are  free  of  discontinuities.  As  an  example, 
the  Hanning  window  produces  a  low  window  error  compared  to  Hamming  and  Taylor 
windows,  which  have  a  discontinuity  at  the  boundaries  of  the  region  of  support.  The 
continuous  Hanning  window  is  a  raised  cosine 


C) 


j  (cos  C  +  l)/2,  c  e 
0,  otherwise, 


(114) 


where  (  is  the  region  of  support.  The  next  two  sub-subsections  analyze  the  window 
error  caused  by  using  the  sum  of  scaled,  shifted,  and  weighted  Hanning  subwinclows  to 
approximate  a  fullband  window  in  frequency  and  a  full  aperture  window  in  azimuth. 


6.2.1  Frequency  Subwindows. 

In  frequency,  the  design  goal  is  to  approximate  a  fullband  Hanning  window  of 
width,  B,  centered  at  fc.  Thus,  the  fullband  Hanning  window  is  succinctly  expressed 


138 


as  Hb  =  H  (J^C  —  fc)-  Due  to  symmetry,  an  odd  number  of  subwindows  produces  a 
smaller  window  error  than  does  an  even  number  of  subwindows.  As  such,  the  number 
of  subwindows,  /,  is  equal  to  2L  +  1,  L  e  N.  Under  these  conditions,  the  error  due  to 
the  frequency  window  approximation  is 


£b  =  " #  (  Jk-/< 


—  ciHb * 


12, 


l=-L 


(115) 


where  HBn  are  the  frequency  subwindows  with  weights,  q.  Here,  the  region  of  sup¬ 
port  for  each  Hanning  subwindow  represents  a  subband  of  the  radar  data.  In  order  to 
obtain  processing  efficiency  in  the  integrated  algorithm,  the  subwindows  in  frequency 
are  scaled  to  exactly  half  of  the  width  of  the  fullband  Hanning  window.  This  choice 
agrees  with  the  fast  CBP  algorithm  presented  in  [12].  Although  additional  process¬ 
ing  efficiency  is  gained  by  narrowing  these  subwindows  even  further,  this  causes  a 
significant  decrease  in  both  classification  and  imaging  accuracy  due  to  the  reciprocal 
relationship  between  bandwidth  and  resolution.  As  the  subbands  narrow,  subimage 
resolution  becomes  coarser,  causing  decreased  localization  in  the  scattering  center 
classifier  and  increased  errors  during  subimage  interpolation. 

Therefore,  the  weighted,  halfband  Hanning  subwindows  in  frequency  are  expressed 
as 

c‘Hm  =  o,H  (Jh  +  (§  -  /„)  .  (116) 

where  the  weights,  q  are  given  by  the  solution  to  the  system  of  equations 


{HB'(-l),HBi(_L))  ■■ 

■■  {Hb>(-L),HB'L) 

C-L 

(Hb,Hb,{_l)) 

(HB'l,  HBi(-l)) 

{Hb>l,HB'B) 

cl 

(Hb,  HB'l) 

(117) 

where  (•,  •)  denotes  the  inner  product  [120].  In  summary,  when  designing  subwindows 
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in  frequency  for  circular  SAR,  the  number  is  odd,  the  spacing  is  given  by  Equation 
(116),  and  the  weights  are  given  by  Equation  (117). 

6.2.2  Azimuth  Subwindows. 

In  azimuth,  the  error  due  to  the  azimuth  window  approximation  is 

j 

ee  =  \\He(e-9c)-J2cjH&j\\2.  (118) 

l=i 

This  error  is  best  controlled,  in  general,  for  an  odd  number,  J,  of  halfway  overlapping, 
uniformly  spaced,  Hanning  subwindows,  so  that 

20 

(119) 

This  condition  determines  the  scaling  of  the  subaperture  windows  so  that 

CjHws  =  CjH  (A<  +  A  ,  (120) 

where  Hqij  are  the  frequency  subwindows  with  weights,  cy  The  weights  and  shifts 
that  reduce  the  error  due  to  the  azimuth  window  approximation  will  depend  upon 
the  type  of  aperture  window  to  be  approximated. 

For  surveillance  SAR  applications  or  ISAR  radar  range  measurements,  360°  aper¬ 
tures  are  common.  In  this  case,  the  full  aperture  window  is  flat,  as  shown  in  the 
bottom  right-hand  corner  of  Figure  47  of  the  intro  to  Chapter  VI.  Fortunately,  the 
use  of  Hanning  subwindows  produces  low  window  approximation  error  in  azimuth 
because  two  equally  weighted  Hanning  windows,  shifted  by  a  difference  of  ir,  sum  to 
unity  in  the  region  of  mutual  support.  This  useful  property  is  a  consequence  of  the 
trigonometric  relationship  (cost)  +  l)/2  +  (cos(C  —  7r)  +  l)/2  =  1,  as  illustrated  in 
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Figure  52.  A  Hanning  window  (left)  has  the  property  that  the  sum  of  equally  weighted, 
halfway  overlapping  Hanning  windows  (right)  is  unity  over  the  region  of  mutual  sup¬ 
port. 


Figure  52.  By  this  property,  a  set  of  properly  shifted  and  scaled  Hanning  subwindows 
of  equal  weighting  will  sum  to  a  flat  window  in  azimuth.  Therefore,  for  a  flat  window, 
Hq  =  1,  Cj  =  1 ,  Vj ,  and  9j  =  in  Equations  (118)  and  (120). 

For  apertures  less  than  360°,  the  full  aperture  window  is  Hq  =  H  (J|£  —  9C ).  In 
this  case,  the  window  approximation  error  is  best  controlled  when  the  windows  are 
uniformly  spaced  by  9j  =  j=^-  —  7 1  —  9C  and  when  the  weights  follow  the  envelope 
of  a  Hanning  window  in  order  to  mimic  the  envelope  of  the  full  aperture  window 
with  Cj  =  H  (j^j  —  7r).  Alternately,  the  weights  can  be  determined  by  adapting  the 
procedure  used  for  subwindows  in  frequency  expressed  in  Equation  (117). 

Finally,  note  that  for  stripmap  Mode  SAR,  the  available  aperture  is  limited  by  the 
antenna  beamwidth,  as  described  in  Section  4.5.  In  fact,  the  available  aperture  is  often 
so  limited  that  only  one  subaperture  can  be  used  per  image.  In  this  degenerate  case, 
J  =  1,  and  no  efficiency  can  be  obtained  by  decomposition  in  the  azimuth  direction. 
However,  subband  windowing  in  frequency  is  still  viable,  and  the  subimages  can  still 
be  weighted  and  summed  to  produce  a  domain  decomposition  image. 

6.2.3  Imaging  Accuracy. 

In  general,  the  imaging  accuracy  improves  as  the  number  of  subwindows  increases. 
For  instance,  frequency  subwindows  for  the  cases  of  /  =  3  and  I  =  5  are  depicted  in 
Figure  53,  where  the  window  error  is  visually  indiscernible  for  the  case  of  /  =  5.  In 
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1  =  3 


1  =  5 


Figure  53.  In  frequency,  a  fullband  Hanning  window  (thick,  dashed  line)  is  approxi¬ 
mated  by  the  sum  (thick,  solid  line)  of  an  odd  number  of  weighted,  shifted,  halfband 
Hanning  windows  (thin,  dotted  lines).  The  weights  for  the  narrow  Hanning  windows 
were  determined  using  Eq.  (117)  for  1  =  3  (left)  and  1  =  5  (right). 


Point  Spread  Function  in  the  Range  Dimension  (dB) 


Figure  54.  Point  spread  functions  in  the  range  dimension  are  oversampled  to  show 
sidelobe  structure.  The  PSF  for  the  case  of  I  =  5  is  visually  similar  to  that  of  the  ideal 
Hanning  function  down  to  approximately  -40  dB.  In  contrast,  the  PSF  for  the  case  of 
I  =  3  is  only  visually  similar  down  to  approximately  -20  dB. 


order  to  determine  how  the  window  approximation  error  will  affect  the  image,  it  is 
helpful  to  calculate  the  point  spread  function  (PSF).  The  PSF  predicts  the  response 
of  a  single,  isotropic,  non-distributed,  non- dispersive  scatterer  located  at  the  zero- 
phase  reference.  In  other  words,  it  illustrates  how  an  ideal  point  scatterer  located 
in  the  center  of  the  image  will  appear.  The  PSF  in  the  range  dimension  is  simply 
the  absolute  value  of  the  inverse  discrete  Fourier  transform  of  Hs{f)  [76].  Figure 
54  compares  the  ideal  PSF  for  a  Hanning  window  to  the  PSFs  for  the  approximated 
Hanning  windows  when  1  =  5  and  1  =  3.  For  the  case  of  /  =  5,  the  PSF  is  visually 
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Window  Error  in  Azimuth 


Figure  55.  The  window  approximation  error  in  azimuth  using  discrete  Hanning  sub¬ 
windows.  For  a  given  aperture,  fewer  subwindows  results  in  more  digital  samples  per 
subwindow. 


similar  down  to  approximately  -40  dB,  but  for  the  case  of  /  =  3,  this  is  only  true  down 
to  approximately  -20  dB.  As  a  result,  the  circular  SAR  image  artifacts  due  to  window 
error  in  frequency  for  the  integrated  algorithm  are  likely  to  be  visually  indiscernible 
for  the  case  of  I  =  5,  but  may  be  noticeable  for  the  case  of  /  =  3.  Section  6.2.3. 1 
presents  empirical  evidence  to  support  these  predictions. 

For  the  case  of  approximating  a  flat  window  in  azimuth,  the  window  error  in 
azimuth  is  zero  for  continuous  Hanning  subwindows,  that  is  e©  =  0.  Consequently, 
for  the  case  of  discrete  Hanning  subwindows,  the  window  error  approaches  zero  as  the 
number  of  samples  for  each  subwindow  approaches  infinity.  Figure  55  illustrates  this 
with  results  produced  using  the  hanning  function  in  Matlab®.  For  a  given  aperture, 
fewer  subwindows  results  in  more  digital  samples  per  subwindow.  Note  also  how  the 
absolute  window  error  decreases  as  the  number  of  samples  for  each  of  the  Hanning 
subwindows  increases.  A  typical  implementation  of  the  integrated  algorithm  features 
dozens  of  subapertures,  each  with  hundreds  of  samples.  In  this  case,  the  per  sample 
window  error  in  azimuth  is  very  low,  and  the  resulting  image  artifacts  are  likely  to 
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Figure  56.  In  azimuth,  a  wide  Hanning  window  (thick,  dashed  line)  can  be  approx¬ 
imated  by  the  sum  (thick,  solid  line)  of  an  odd  number  of  weighted,  shifted  narrow 
Hanning  windows  (thin,  dotted  lines).  The  approximation  error  is  reduced  as  the 
number  of  narrow  windows,  J,  increases. 

be  visually  indiscernible.  Section  6.2.3. 1  presents  empirical  evidence  to  support  this 
prediction. 

In  contrast,  for  the  case  of  approximating  a  Hanning  window  in  azimuth,  the 
imaging  accuracy  will  improve  as  the  number  of  subwindows  increases,  as  indicated 
in  Figure  56.  The  PSFs  for  this  case  are  shown  in  Figure  57.  Note,  that  there  must 
be  at  least  J  =  7  subwindows  in  order  to  keep  the  sidelobes  in  azimuth  below  -30  dB. 
Recall  that  the  wide-angle  approximation  in  the  multi-peak  model  requires  0'  >  10°. 
Also  note  that  according  to  Equation  (119),  a  full  aperture  of  40°  is  required  in 
order  to  create  J  =  7  subapertures  of  width  0'  =  10°  each.  This  suggests  that  the 
integrated  algorithm  requires  a  very  wide  aperture,  0  >  40°,  in  order  to  reduce  image 
artifacts  to  a  visually  indiscernible  level. 

Last,  note  that  the  subwindow  design  principles  presented  earlier  are  not  unique  to 
Hanning  windows  and  are  easily  adapted  for  use  with  other  types  of  windows.  How¬ 
ever,  it’s  worth  restating  that  windows  with  no  discontinuities,  such  as  the  Hanning 
window,  produce  a  lower  window  error  than  other  popular  windows  which  feature 
discontinuities,  such  as  Hamming  and  Taylor  windows. 
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Figure  57.  The  PSFs  in  the  cross-range  dimension  are  oversampled  to  show  sidelobe 
structure.  The  sidelobe  levels  in  cross-range  are  reduced  as  the  number  of  subaperture 
windows,  J,  increases. 


6.2.3. 1  Empirical  Studies  of  Imaging  Accuracy. 

This  section  compares  the  imaging  accuracy  of  the  integrated  algorithm  to  that 
of  conventional  CBP  for  a  surveillance  SAR  application.  As  previously  discussed,  an 
image  produced  by  the  integrated  algorithm  contains  image  artifacts  due  to  errors  in 
interpolation  and  window  approximation.  Additional  analysis  revealed  that  the  errors 
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Figure  58.  SAR  images  of  a  parking  lot  scene  (bottom  row)  with  call-out  boxes  fea¬ 
turing  top-hat  and  dihedral  calibration  targets  (top  row).  The  images  produced  by 
conventional  CBP  (left  column)  serve  as  a  baseline  for  evaluating  the  imaging  accuracy 
of  the  integrated  algorithm.  The  errors  introduced  by  the  integrated  algorithm  for 
the  case  of  I  =  5  (middle  column)  are  imperceptible.  However,  the  errors  introduced 
for  the  case  of  I  =  3  (right  column)  cause  a  visible  increase  in  sidelobe  energy,  which 
causes  the  appearance  of  a  slightly  stronger  and  broader  halo  around  objects  with  high 
scattering  intensity. 


were  quite  small  and  likely  to  be  visually  indiscernible  for  the  case  of  a  flat  window 
in  azimuth  and  1  =  5  subbands.  However,  for  the  case  of  /  =  3,  the  error  due  to 
window  approximation  in  frequency  is  appreciable  and  is  likely  to  be  noticeable  due  to 
increased  PSF  sidelobes.  The  following  examples  verify  these  conclusions  empirically. 

Figure  58  shows  SAR  images  of  the  Gotcha  data  described  in  Table  15  of  Section 
5.4.4.  Note  that  the  aperture  is  360°  so  that  a  flat  aperture  window  is  approximated 
in  azimuth.  For  brevity,  only  the  data  from  the  horizontal  transmit,  horizontal  receive 
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(HH)  channel  is  shown.  Also,  in  order  to  emphasize  differences  between  imaging 
methods,  a  coherent  subaperture  sum  with  2x  oversampling  is  used.  The  SAR  images 
are  reconstructed  using  conventional  CBP,  the  integrated  algorithm  with  I  =  5, 
and  the  integrated  algorithm  with  1  =  3.  The  call-out  boxes  in  the  top  row  reveal 
additional  detail  for  the  top  hat  and  one  trihedral  taken  from  the  full-scene  images  in 
the  bottom  row.  By  using  the  conventional  CBP  images  (left  column)  as  a  baseline, 
one  can  visually  confirm  the  impact  of  the  image  artifacts  introduced  by  the  integrated 
algorithm.  For  the  case  of  /  =  5  (middle  column),  the  artifacts  are  imperceptible, 
and  the  images  in  the  left  and  middle  columns  are  virtually  identical.  In  contrast,  for 
the  case  of  /  =  3  (right  column),  some  image  artifacts  are  noticeable.  The  increase 
in  PSF  sidelobe  energy  causes  a  slightly  stronger  and  broader  halo  to  appear  around 
objects  with  high  scattering  intensity.  As  a  result,  some  of  the  pixel  values  are  visibly 
different  when  compared  to  conventional  CBP,  especially  those  in  the  center  and 
along  the  edge  of  the  top  hat  and  those  surrounding  the  trihedral.  Similar  results 
(not  shown)  were  obtained  for  varying  amounts  of  oversampling  in  the  integrated 
algorithm. 

Results  in  Figure  58  are  valuable  for  analyzing  the  usefulness  of  the  integrated 
algorithm  and  for  gaining  practical  insight  into  the  engineering  trade  space  between 
imaging  accuracy  and  processing  efficiency.  For  the  case  of  /  =  5,  the  imaging  accu¬ 
racy  is  very  good  and  any  processing  gain  makes  it  a  superior  choice  to  conventional 
CBP.  However,  for  the  case  of  /  =  3,  the  SAR  system  engineer  must  evaluate  the 
noticeable  loss  in  imaging  accuracy  against  the  benefits  of  associated  processing  gain. 
The  cases  for  /  >  5  are  not  presented  because  the  processing  cost  increases  linearly 
with  no  appreciable  improvement  in  imaging  accuracy.  The  next  section  analyzes  the 
processing  cost  of  the  integrated  algorithm  and  demonstrates  its  potential  efficiency. 
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6.3  Efficiency  Study 


Scattering  center  classification  by  phase  history  decomposition  is  operationally 
efficient  when  compared  to  existing  image  segmentation  methods.  This  is  because 
image  segmentation  requires  human  supervision,  as  discussed  in  Sections  1.2  and 
3.1.  Although  the  expected  gains  in  operational  efficiency  are  likely  to  be  substantial, 
these  are  difficult  to  analyze  mathematically.  In  contrast,  the  computational  efficiency 
gains  afforded  by  the  integrated  algorithm  are  more  readily  analyzed. 

The  integrated  algorithm  allows  gains  in  computational  efficiency  with  a  control¬ 
lable  reduction  in  imaging  accuracy.  Use  of  this  trade-space  is  desirable  for  highly  re¬ 
source  constrained  scenarios,  such  as  SAR  surveillance  applications.  Imaging  accuracy 
was  discussed  in  the  preceding  sections,  while  computational  efficiency  is  discussed 
in  this  section.  The  section  begins  by  introducing  the  concept  of  SAR  surveillance 
and  the  processing  constraints  that  limit  maximum  coverage  area.  Then,  it  presents 
analytical  analysis  of  the  expected  efficiency  of  the  integrated  algorithm  assuming 
a  convolution  backprojection  imaging  algorithm.  Finally,  it  uses  coverage  area  as 
a  metric,  to  present  empirical  studies  of  computational  efficiency  for  the  integrated 
algorithm. 

6.3.1  SAR  Surveillance. 

SAR  systems  are  useful  for  wide-area  surveillance,  provided  that  these  systems 
can  produce  images  in  near  real-time  [118,  18,  99,  25].  Example  applications  include 
all-weather,  day  and  night  monitoring  of  borders,  roads,  or  cities  to  protect  popula¬ 
tions  from  illicit  activities.  As  a  result  of  the  growing  demand  for  SAR  surveillance 
imagery,  developers  of  airborne  SAR  systems  have  begun  to  incorporate  a  extensive 
circular  SAR  modes  into  their  designs  [23,  22,  35,  151,  89].  Surveillance  via  circu¬ 
lar  SAR  imaging  is  depicted  in  Figure  2  of  Section  2.1.1.  The  operating  parameters 
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and  capabilities  of  specific  surveillance  SAR  systems  vary,  but  typical  performance 
parameters  include  a  resolution  of  0.15-m  to  0.5-m  and  a  coverage  area  of  0.5-km2 
to  1-km2.  Current  efforts  seek  to  increase  the  available  coverage  area  for  a  given  res¬ 
olution  under  time  and  resource  constraints.  The  maximum  coverage  area  depends 
on  the  complexity  of  the  imaging  algorithm,  the  speed  of  the  image  processor,  and 
the  computational  burden  for  optional  post-processing  of  the  SAR  data,  and  serves 
as  the  performance  metric  for  the  efficiency  studies  in  Section  6.3.3. 

First,  coverage  area  depends  upon  the  complexity  of  the  imaging  algorithm.  The 
imaging  algorithms  suitable  for  circular  SAR  are  CBP  and  PFA  [76].  In  theory, 
both  algorithms  produce  equivalent  images  with  controllable  error  [82],  Furthermore, 
for  an  image  consisting  of  N  x  N  pixels,  each  has  a  computational  complexity  of 
0(N2  log  N)  when  fast  Fourier  transforms  and  fast  CBP  techniques  are  considered 
[103,  12].  Last,  both  algorithms  are  scalable  to  allow  accelerated  processing  through 
distributed  and  parallel  architectures  [82,  68].  However,  despite  these  and  other  gen¬ 
eral  similarities,  CBP  has  one  significant  advantage  over  PFA:  it  allows  for  greater 
flexibility  in  selecting  the  image  pixel  locations  and  spacing,  which  can  be  used  to  fo¬ 
cus  the  image  to  a  previously-obtained  digital  elevation  map  [76,  40].  For  this  reason, 
the  CBP  algorithm  is  generally  preferred  in  surveillance  SAR  applications,  where  the 
desired  imaging  plane  often  varies  significantly  with  changes  in  aspect  angle.  The 
computational  complexity  varies  from  0(N 3)  for  conventional  CBP  to  0(N2  log  N) 
for  fast  CBP  techniques,  depending  upon  the  implementation.  As  imaging  algorithms 
are  mature,  significant  increases  in  real-time  (or  near  real-time)  coverage  area  will  be 
achieved  by  increased  processing  capability  and  efficiency. 

Second,  coverage  area  depends  upon  the  speed  of  the  image  processor.  Advances  in 
parallel  processors  have  lowered  cost  and  increased  performance  of  near  real-time  SAR 
imaging  systems  [65,  68,  41].  However,  the  coverage  area  of  current  systems  is  still 
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quite  limited.  For  example,  Gorham  et  al.  reported  on  the  processing  time  required 
to  produce  circular  SAR  images  using  a  specialized  processor  suitable  for  surveillance 
aircraft  [65] .  The  processor  produces  a  2048  x  2048-pixel  circular  SAR  image  for  one 
polarization  channel  in  approximately  23  seconds.  This  provides  a  0.25-km2  coverage 
area  for  a  single  polarization  channel  at  0.25-m  resolution.  Assuming  systematic 
scaling,  the  processor  produces  a  fully-polarized  circular  SAR  image  with  similar 
resolution  over  a  4-km2  coverage  area  in  approximately  45  minutes.  For  comparison, 
a  separate  implementation  using  four  graphics  processing  units  produces  the  same 
result  in  approximately  90  minutes  [68].  Because  a  typical  airborne  radar  system 
circumnavigates  its  entire  flight  path  in  under  45  minutes,  a  4-km2  coverage  area 
exceeds  what  these  solutions  can  provide  in  real-time.  As  processing  advances  are 
incremental,  even  with  fast  imaging  algorithms,  increases  in  real-time  coverage  area 
are  expected  to  remain  incremental. 

Third,  the  information  within  the  coverage  area  depends  upon  the  computational 
burden  for  optional  post-processing  of  the  SAR  data.  Upon  receipt  of  a  SAR  image, 
an  analyst  must  subsequently  scrutinize  the  image,  but  often  requires  automatic  anal¬ 
ysis  tools  to  meet  time  constraints.  The  time  required  for  image  analysis  is  difficult 
to  parameterize.  However,  it  is  generally  accepted  that  people  can  better  prioritize 
surveillance  resources  and  improve  image  analysis  with  the  aid  of  additional  SAR 
data  processing  [116].  Proposed  methods  range  from  simple  scatterer  filters  or  clas¬ 
sifiers,  which  help  identify  regions  containing  objects  of  interest,  to  automatic  target 
classifiers,  which  provide  specific  target  information  [116].  While  the  advantages  and 
limitations  of  such  classifiers  vary,  almost  all  of  them  execute  significant  amounts  of 
signal  processing  after  a  SAR  image  is  formed  [116].  For  instance,  some  advanced 
methods  employ  iterative  model-matching  techniques,  which  require  repetitive  imag¬ 
ing  of  the  scene  [101,  122],  In  addition,  many  signal  estimation  methods  require 
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evaluation  of  a  full-rank  correlation  matrix,  which  drives  computational  complexity 
as  high  as  0(A6)  [39].  Unfortunately,  a  dramatic  increase  in  processing  costs  occurs 
when  such  classifiers  are  appended  to  a  SAR  imaging  process.  While  researchers  of¬ 
ten  propose  efficient  implementations  for  classifiers  with  high  orders  of  computational 
complexity,  these  are  tailored  to  the  particular  classifier,  and  research  to  survey  the 
processing  costs  of  all  such  classifiers  is  nascent.  As  a  result,  current  evaluations  of 
the  near-real  time  capability  of  circular  SAR  systems  do  not  usually  consider  op¬ 
tional  post-processing  of  the  SAR  data  as  part  of  the  equation.  Therefore,  advances 
in  optional  post-processing  of  the  SAR  data  alone  have  not  translated  into  increased 
coverage  area  for  circular  SAR  systems,  to  date.  However,  the  integrated  algorithm 
has  the  potential  to  reduce  or  even  eliminate  the  post-processing  costs  of  some  clas¬ 
sifiers  and  produce  a  significant  increase  in  coverage  area  for  circular  SAR  systems, 
as  discussed  in  Sections  1.2  and  3.1. 

6.3.2  Analysis  of  Computational  Cost. 

Because  domain  decomposition  imaging  has  a  lower  order  of  computational  com¬ 
plexity  than  conventional  CBP  imaging,  the  integrated  algorithm  has  the  potential 
to  produce  a  SAR  image  annotated  with  scatterer  classification  in  less  time  than  it 
takes  to  produce  the  same  image  by  conventional  CBP,  without  scatterer  classifica¬ 
tion.  This  section  analyzes  the  computational  costs  of  both  the  integrated  algorithm 
and  conventional  CBP  to  determine  under  what  conditions  this  efficiency  is  achieved. 

A  first-order  analysis  of  computational  cost  scales  linearly  with  an  increase  in  the 
number  of  available  polarization  channels.  Therefore,  the  following  analysis  is  re¬ 
stricted  to  a  single  polarization  channel.  Also,  to  facilitate  discussion,  primed  math¬ 
ematical  constants  refer  to  the  integrated  algorithm,  while  unprimed  constants  refer 
to  conventional  CBP.  Note  that  the  equation  for  conventional  CBP  imaging  using 
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overlapping  subapertures  is  given  by  Equation  (4)  of  Section  2. 1.2.4,  whereas  the  in¬ 
tegrated  algorithm  uses  phase  history  decomposition  as  expressed  in  Equation  (3)  of 
Section  2. 1.2.2.  For  both,  the  imaging  operator  has  a  computational  complexity  of 
0(N 3)  for  an  N  x  N  image,  while  other  processes  have  a  computational  complexity 
of  0(N2)  [76], 

The  computational  cost  for  CSAR  image  formation  by  conventional  CBP  is 

Cost  =  J(bN3  +  sN2),  (121) 

where  b  is  an  unspecified  coefficient  for  the  imaging  operator  and  s  is  an  unspecified 
coefficient  for  other  processes.  These  0(N2)  processes  include  subaperture  image 
summation  and  various  data  conditioning  operations.  Examples  of  data  conditioning 
operations  are  two-dimensional  (2D)  windowing  of  the  subaperture  data  and  compu¬ 
tation  of  pixel  intensities  for  use  in  displaying  an  image. 

In  contrast,  the  integrated  algorithm  creates  potential  savings  by  producing  subim¬ 
ages  that  are  |  x  f  in  size  while  adding  an  additional  2D  interpolation  process  and 
a  feature  extraction  process  both  of  which  are  0(N2).  The  number  of  subimages  is 
/  x  J,  and  the  computational  cost  to  produce  an  annotated  CSAR  image  using  the 
integrated  algorithm  is 

Cost'  =  If  (V  (f  )3  +  s'  (f  )2  +  a'  (f  )2)  ,  (122) 

where  a'  is  an  unspecified  coefficient  which  accounts  for  the  additional  2D  inter¬ 
polation  and  feature  extraction  processes,  and  where  b'  and  s'  are  the  fast  CBP 
counterparts  to  b  and  s. 

Assuming  overlapping  subapertures  and  subimage  dimensions  with  equal  resolu¬ 
tion,  the  widths  of  the  subapertures  for  conventional  CBP  are  twice  the  width  of 
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the  subapertures  for  the  integrated  algorithm.  This  occurs  because  the  integrated 
algorithm  recovers  subimages  from  halfband  data,  where  B'  —  As  a  result,  the 
integrated  algorithm  requires  twice  as  many  subaperture  images  as  does  conventional 
CBP,  and  J'  =  2 J.  It  can  be  shown  under  the  preceding  assumptions,  that  the 
relationship 

N  >  2I(s'  +  a)  -4s  Cost'  <  Cost  (123) 

46-/6'  ~ 

holds.  When  I  >  3,  the  numerator  is  likely  to  be  positive.  Therefore,  as  long  as 
46  >  lb',  there  exists  some  threshold,  N  >  Nthresh ,  at  which  the  cost  of  integrated 
algorithm  is  less  than  conventional  CBP.  The  conditions  which  cause  the  denominator 
of  (123)  is  positive,  warrant  further  scrutiny.  Intuitively,  it  seems  safe  to  assume  that 
6  >  6'  because  the  subimages  are  of  a  coarser  resolution  for  the  integrated  algorithm. 
Based  on  this  assumption,  Nthresh,  exists  for  the  case  of  /  =  3.  However,  it  is  unclear 
whether  or  not  Nthresh  exists  for  the  case  of  /  >5.  In  general,  Nthresh  must  be 
determined  empirically  because  a',  6,  and  6'  are  unspecified.  For  instance,  6  is  affected 
by  the  computational  complexity  of  the  interpolation  in  the  imaging  operator,  £?{•}, 
while  a'  is  affected  by  the  computational  complexity  of  the  2D  interpolation  in  the 
interpolation  operator,  Z{-}. 

In  summary,  first-order  comparison  of  the  computational  cost  of  the  integrated 
algorithm  to  that  of  conventional  CBP  showed  that  for  large,  high-resolution  images 
with  large  N,  the  integrated  algorithm  has  a  lower  processing  cost  than  conventional 
CBP  for  the  case  of  /  =  3,  but  not  necessarily  for  the  case  of  /  >  5.  In  addi¬ 
tion,  the  integrated  algorithm  produces  scatterer  classification  information,  whereas 
conventional  CBP  does  not.  The  next  section  presents  empirical  studies  of  the  net 
computational  savings  and  imaging  accuracy  of  the  integrated  algorithm. 
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Figure  59.  Processing  times  required  to  produce  a  circular  SAR  image.  For  large  TV, 
conventional  CBP  has  a  higher  processing  cost  than  the  integrated  algorithm  for  the 
case  of  I  =  3. 


6.3.3  Empirical  Studies  of  Processing  Efficiency. 

This  section  evaluates  the  processing  efficiency  of  the  integrated  algorithm,  as 
compared  to  conventional  CBP.  The  discussion  in  Section  6.3.2  revealed  that  there 
is  a  theoretical  image  size,  N  >  Nthresh}  above  which  the  integrated  algorithm  has  a 
lower  processing  cost  than  conventional  CBP.  Additional  analysis  revealed  that  this 
threshold  exists  for  the  case  when  I  =  3,  but  it  is  unclear  if  it  exists  for  the  case  when 
1  =  5.  The  following  examples  verify  these  conclusions  empirically. 

For  near  real-time  surveillance,  the  time  required  to  produce  an  N  x  N  image  is  of 
interest.  The  following  empirical  studies  use  the  Gotcha  data  set  described  in  Table 
15  of  Section  5.4.4  to  simulate  a  SAR  surveillance  scene.  To  simulate  different  values 
of  N,  the  number  of  samples  in  the  Gotcha  data  set  is  adjusted  by  an  appropriate 
amount  of  decimation  or  interpolation.  Based  on  this  procedure,  Figure  59  reports 
the  processing  times  for  a  series  of  simulations  of  varying  image  sizes.  The  longer 
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run-times  were  approximated  by  processing  only  a  portion  of  the  data  set.  For  the 
case  of  /  =  3,  the  integrated  algorithm  has  a  lower  processing  time  than  conven¬ 
tional  CBP  for  N  >  200.  In  contrast,  for  I  =  5,  the  integrated  algorithm  always 
has  a  higher  processing  time  than  conventional  CBP.  For  large  N,  all  three  curves 
display  a  computational  complexity  of  0(N3),  as  expected.  These  simulations  were 
conducted  using  a  64-bit  version  of  Matlab®  on  a  quad-core  Intel®  Xeon®  5160, 
3-GHz  processor. 

A  useful  metric  for  evaluating  the  integrated  algorithm  in  a  time  constrained  sce¬ 
nario  is  the  amount  of  coverage  area  either  gained  or  lost,  as  compared  to  conventional 
CBP.  Figure  60  presents  this  comparison  for  the  three  cases  shown  in  Figure  59.  To 
construct  Figure  60,  it  is  noted  that  the  curves  in  Figure  59  were  interpolated  by  a 
piecewise  cubic  hermite  interpolating  polynomial  provided  in  the  pchip  function  of 
Matlab®.  Based  on  this  procedure,  Figure  60  shows  the  estimated  coverage  area  for 
the  integrated  algorithm  as  a  percent  of  the  coverage  area  for  conventional  CBP.  The 
coverage  area  for  the  integrated  algorithm  is  approximately  20-percent  greater  than 
conventional  CBP  for  the  case  of  /  =  3  and  about  10-percent  less  than  conventional 
CBP  for  the  case  of  /  =  5. 

The  best  results  for  the  integrated  algorithm  occur  for  N  «  565,  and  so  the 
performance  of  the  algorithm  near  this  operating  point  deserves  additional  evaluation. 
Of  particular  interest  from  Equation  (122)  is  the  condition  under  which  4 b  >  5 b' , 
where  the  processing  cost  of  conventional  CBP  is  greater  than  the  processing  cost  of 
the  integrated  algorithm  (shown  for  the  case  of  /  =  5).  The  relationship  between 
b  and  b'  is  primarily  affected  by  the  processing  cost  of  the  one-dimensional  (ID) 
interpolation  in  the  imaging  operator.  For  each  backprojection,  the  image  pixels 
accumulate  interpolated  samples  of  the  ID  range-compressed,  filtered  radar  data  [76]. 
In  order  to  obtain  acceptable  imaging  accuracy  and  efficiency,  the  imaging  operator 
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Figure  60.  Estimated  circular  SAR  coverage  areas  for  the  integrated  algorithm  as  a 
percent  of  the  coverage  area  for  conventional  CBP.  For  large  N,  where  processing  times 
are  large,  the  integrated  algorithm  provides  approximately  20-percent  more  coverage 
than  conventional  CBP  when  1  =  3,  and  it  provides  approximately  10-percent  less 
coverage  when  1  =  5. 

usually  performs  ID  interpolation  by  a  two-step  process  of  ideal  sine  interpolation 
followed  by  a  much  less  computationally  complex  linear  interpolation  [64] .  In  practice, 
ideal  sine  interpolation  is  accomplished  by  zero-padding  the  radar  data  in  the  spectral 
domain  before  performing  an  IFFT  to  obtain  an  oversampled  signal  in  the  spatial 
domain.  Typically,  the  amount  of  oversampling  is  on  the  order  of  10  x,  or  more 
[64,  135].  For  instance,  approximately  9x  oversampling  was  used  in  these  studies. 

Additional  experiments  revealed  that  simply  reducing  oversampling  in  both  con¬ 
ventional  CBP  and  the  integrated  algorithm  to  4x  caused  the  computational  cost  of 
traditional  CBP  to  become  nearly  equivalent  to  the  processing  cost  of  the  integrated 
algorithm  for  the  case  of  /  =  5  and  N  =  565.  However,  reducing  oversampling  to 
such  a  low  amount  produced  image  artifacts  that  are  usually  undesirable.  Therefore, 
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these  studies  suggest  that  for  typical  SAR  surveillance  implementations,  the  process¬ 
ing  cost  for  conventional  CBP  will  be  less  than  the  processing  cost  for  the  integrated 
algorithm  for  the  case  of  /  =  5. 

Finally,  note  that  these  efficiency  studies  assumed  that  the  integrated  algorithm 
required  twice  as  many  subaperture  images  as  conventional  CBP,  where  J'  =  2  J  in 
Equation  (123).  If  this  requirement  is  relaxed,  so  that  the  subaperture  width  for 
the  integrated  algorithm  is  equal  to  that  for  conventional  CBP,  then  computational 
efficiency  will  improve.  However,  the  efficiency  is  not  likely  to  double  when  J'  =  J 
because  larger  subapertures  may  demand  a  higher  order  interpolator,  as  discussed 
previously  in  Section  6.1.1.  However,  when  J'  =  J,  it  is  likely  that  there  exists  a 
threshold,  Nthresh ,  whereby  the  coverage  area  is  larger  for  both  cases  of  /  =  3  and 
I  =  5  when  compared  to  the  coverage  area  for  conventional  CBP. 
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VII.  Conclusion 


In  conclusion,  this  dissertation  presented  a  new  theory  for  SAR  scatterer  classifi¬ 
cation  by  phase  history  decomposition.  More  specifically,  the  hypotheses  presented  in 
Section  3.2  have  been  proven.  First,  canonical  scatterers  can  be  located  and  classified 
in  subimages.  This  is  accomplished  using  the  multi-peak  model  and  SPLIT  algorithm. 
Second,  the  process  is  efficient.  Operational  efficiency  is  improved  by  eliminating  the 
need  for  human  supervision  associated  with  existing  scatterer  classification  methods 
based  on  image  segmentation.  Furthermore,  computational  efficiency  is  improved 
through  use  of  the  integrated  algorithm. 

The  benefits  of  the  new  theory  are  significant.  First,  it  provides  a  solid  theoret¬ 
ical  basis  by  which  to  predict  imaging  and  classification  performance.  Second,  it  is 
flexible,  with  a  controllable  efficiency  versus  accuracy  trade-space  that  can  be  used 
to  optimize  desired  performance  for  resource  constrained  scenarios.  Third,  the  new 
theory  provides  a  general  tool  for  queuing  analysts  or  precision  algorithms  for  high 
precision,  post-processing  of  SAR  images.  Finally,  note  that  this  research  has  resulted 
in  the  publication  of  three  conference  papers,  the  submission  of  two  journal  articles, 
and  the  production  of  computer  programs  for  use  with  real  SAR  data. 

7.1  Notable  Limitations 

While  there  are  many  significant  benefits  of  the  new  theory  for  SAR  scatterer 
classification  by  phase  history  decomposition,  there  are  some  notable  limitations  to 
keep  in  mind  as  well.  The  notable  limitations  of  the  multi-peak  model  are  listed 
in  decreasing  order  of  importance.  First,  model  accuracy  is  limited  by  a  wide-angle 
approximation  that  restricts  apertures  or  subapertures  to  be  greater  than  or  equal 
to  10°.  Fortunately,  wide-angle  SAR  imaging  systems  do  exist  and  their  continued 
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development  makes  wide-angle  SAR  applications  likely  to  be  common  in  the  future. 
Second,  for  stripmap  SAR  applications,  the  multi-peak  model  is  restricted  for  use 
with  the  Omega-K  algorithm.  Fortunately,  the  Omega-K  algorithm  is  mature  and 
even  commonly  used  for  P-band  stripmap  SAR  systems.  Third,  also  for  stripmap 
SAR,  the  multi-peak  model  is  accurate  for  scatterers  with  effective  lengths  greater 
than  or  equal  to  1.6  times  the  physical  length  of  the  antenna.  Fourth,  a  small-angle 
approximation  restricts  the  multi-peak  model  to  apertures  or  subapertures  less  than 
20°  for  rectangularly  formatted  phase  histories.  Although  rectangularly  formatted 
phase  histories  are  very  common  for  spotlight  mode  SAR,  this  restriction  is  already 
often  observed,  in  practice.  Fifth,  the  multi-peak  model  is  limited  by  a  high-frequency 
approximation  for  perfect  electrical  conductors  in  the  far  held.  In  particular,  this 
requires  that  the  effective  length  of  scatterer  to  be  ten  wavelengths  long  or  longer. 
Fortunately,  these  conditions  are  met  for  many  objects  of  interest  in  typical  SAR 
applications.  Furthermore,  these  restrictions  are  common  among  existing  scatterer 
models  and  classification  methods,  and  therefore,  are  not  unique  to  the  multi-peak 
model.  Sixth,  the  multi-peak  model  requires  the  use  of  a  tapered  window  in  azimuth. 
Fortunately,  windowing  of  the  phase  history  is  already  a  common  practice  in  SAR 
imaging. 

The  notable  limitations  of  the  SPLIT  algorithm  are  listed  in  decreasing  order  of 
importance.  First,  classification  accuracy  is  limited  by  subimage  resolution,  which 
is  half  of  the  system  resolution.  The  SPLIT  algorithm  obtains  the  required  spectral 
information  by  reducing  spatial  resolution.  As  expected,  such  trade-offs  are  unavoid¬ 
able  in  SAR  image  processing  clue  to  the  principle  of  time-frequency  reciprocity. 
Second,  the  SPLIT  algorithm  is  generally  inaccurate  for  fractional  bandwidths  be¬ 
low  10-percent.  The  SIR  due  to  neighboring  scatterers  must  be  well-controlled  to 
ensure  good  classification  accuracy.  Experimental  results  indicate  that  a  fractional 
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bandwidth  of  at  least  10-percent  is  needed  to  reduce  the  SIR  for  typical  man-made 
structures.  Third,  the  least  squares  classifier  is  not  optimized  for  clutter  or  interfer¬ 
ence  limited  scenarios.  It  is  well-suited  for  classifying  canonical  scatterers  on  simple 
targets  in  free-space,  but  when  interference  from  neighboring  strong  scatterers  or 
clutter  are  a  concern,  the  performance  of  the  least  squares  classifier  is  expected  to 
be  suboptimal.  Fourth,  the  SPLIT  algorithm  requires  a  normalized  phase  history. 
Fortunately,  the  necessary  normalization  factors  are  typically  available  for  most  SAR 
systems.  Fifth,  the  SPLIT  algorithm  assumes  a  stationary  scattering  center.  As 
expected,  some  objects,  such  as  resonant  cavities,  and  moving  vehicles  or  windmill 
blades  appear  unfocused  in  the  SAR  image.  Fortunately,  these  objects  are  typically 
rejected  by  SPLIT  algorithm  using  the  stable  peak  criterion  and  are  not  expected  to 
cause  undue  confusion  in  the  classifier. 

There  are  three  notable  limitations  of  the  integrated  algorithm.  Fortunately,  if 
these  limitations  prove  unacceptable  in  practice,  then  scatterer  classification  can  be 
accomplished  separately  and  results  can  be  overlaid  on  a  SAR  image  reconstructed 
using  a  conventional  imaging  process.  First,  rapid  phase  variations  in  the  subimages 
require  either  a  high  order  interpolator  or  restricted  subdomain  in  order  to  reduce 
image  artifacts  to  imperceptible  levels.  Second,  the  integrated  algorithm  is  limited  to 
the  use  of  Hanning  windows  to  limit  discontinuities  which  increase  window  approxi¬ 
mation  error.  Third,  the  window  approximation  error  can  introduce  noticeable  image 
artifacts  for  apertures  under  40°  and  a  number  of  subbands  less  than  five. 

7.2  Future  Work 

This  dissertation  creates  a  solid  foundation  for  future  research  in  the  following 
areas.  First,  recent  trends  in  legal  restrictions  on  use  of  the  electromagnetic  spectrum 
have  spurred  interest  in  SAR  collections  over  discontinuous  frequencies.  Fortunately, 
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the  new  theory  in  this  dissertation  can  be  extended  for  use  with  discontinuous  phase 
histories,  as  needed.  Second,  classification  accuracy  can  be  improved  beyond  that  of  a 
least  squares  classifier  by  use  of  standard  signal  detection  and  estimation  techniques. 
For  example,  the  use  of  generalized  likelihood  ratio  tests  for  specific  scenarios  could  be 
of  interest.  Alternately,  adaptive  filters  for  clutter  cancellation  can  be  incorporated 
into  the  algorithm.  Third,  the  theory  can  be  extended  beyond  monostatic,  2D  SAR 
to  include  3D  SAR  and  bi-static  SAR.  Fourth,  domain  decimation  imaging  featuring 
a  mosaic  of  subimages  can  provide  computational  efficiency  equivalent  to  domain 
decomposition  imaging.  These  techniques  can  be  blended  rather  easily,  but  the  effect 
on  classification  accuracy  is  uncertain  and  would  need  to  be  researched.  Finally, 
there  are  potentially  additional  uses  for  related  versions  of  the  multi-peak  model  and 
SPLIT  algorithm  in  other  topics  and  areas  of  research,  such  as  the  CLEAN  algorithm 
modified  by  the  stable  peak  criterion  [53,  100,  141]. 
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Appendix  A.  AFIT  Indoor  RCS  Measurement  Range 


The  AFIT  indoor  RCS  measurement  range  consists  of  a  large  anechoic  chamber, 
target  support  system,  and  a  radar  system  featuring  a  Lintek  gating  box.  Specific 
information  regarding  these  systems  can  be  found  in  the  course  notes  for  EENG  627  - 
RCS  Analysis,  Measurement,  and  Reduction  [34],  Diagrams  of  the  range  are  provided 
in  Figure  61. 
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Side  View 
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Figure  61.  AFIT  indoor  RCS  measurement  range  [104]. 
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